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ABSTRACT 


A  numerical  investigation  of  unsteady  wind  tunnel  and  ground  interference  effects  is  carried 
out  in  the  time  domain  to  study  the  transonic  flutter  characteristics  of  the  NLR  7301  section  inside 
a  wind  tunnel  and  the  thrust  generation  characteristics  of  a  NACA  0014  airfoil  plunging  near  a 
ground  plane.  A  parallelized,  multi-block,  deforming  grid,  unsteady  flow-solver  is  coupled  with  a 
two-degree-of-freedom  structural  model. 

For  the  transonic  flutter  problem,  two  types  of  porous-wall  boundary-conditions  arc  imple¬ 
mented  and  tested  for  the  boundaries  representing  the  tunnel  walls.  The  type  of  porous  boundary 
condition  is  found  to  influence  significantly  both  steady  and  unsteady  solutions.  Results  show  that 
the  free-flight  flutter  behavior  may  differ  significantly  from  the  behavior  found  in  a  porous  wind 
tunnel  because  of  the  strong  dependence  on  the  tunnel  porosity  parameter  and  the  proximity  of  the 
walls. 

An  analysis  of  the  trailing  edge  boundary  condition  is  performed  for  the  airfoil  in  ground 
effect.  The  computations  show  that  this  boundary  condition  influences  the  solution  only  when 
non-linearities  arc  present  in  the  flow-field,  although  parameters  averaged  through  a  cycle  of  os¬ 
cillation  arc  not  affected  significantly.  The  same  behavior  is  observed  for  the  influence  of  the  tur¬ 
bulence  model  on  the  fully-turbulent,  unsteady  computations.  However,  the  best  agreement  with 
low  Reynolds  number,  experimental  data  is  obtained  when  the  flow  is  assumed  laminar  and  no 
turbulence  model  is  applied. 
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I.  INTRODUCTION 


The  study  of  unsteady  flow  problems  has  become  very  important  in  the  past  few  decades. 
Among  many  other  problems,  one  can  mention  thrust  and  lift  generation  by  means  of  oscillating 
wings,  flutter,  dynamic  stall,  blade-row  interactions  in  axial  compressors,  and  blade-vortex  interac¬ 
tions  on  helicopter  blades  as  examples  of  important  unsteady  aerodynamic  problems.  The  problems 
of  transonic  flutter  of  an  airfoil  inside  a  wind  tunnel  and  flapping  wing  propulsion  of  opposed- 
plunging  airfoils  arc  the  primary  interest  of  this  work. 

The  main  purpose  of  this  introductory  chapter  is  to  give  the  reader  an  overview  of  the  latest 
achievements  related  to  these  two  problems  and  to  summarize  the  contributions  of  the  present  work 
to  the  numerical  computation  of  unsteady  flows  past  a  fluttering  airfoil  inside  a  wind  tunnel  and 
over  a  flapping  airfoil  in  ground  effect.  Consequently,  this  chapter  is  divided  into  specific  sections 
related  to  these  two  problems.  It  also  includes  a  section  about  some  of  the  details  of  an  unsteady 
aerodynamics  solver  common  to  both  problems. 

A.  TRANSONIC  FLUTTER  STUDIES 

Historically,  great  efforts  have  been  devoted  to  experimental  and  theoretical  investigations 
of  the  transonic  flutter  characteristics  of  airfoils  because  the  transonic  dip  associated  with  the  flutter 
speed  of  typical  aircraft  wings  poses  a  serious  problem  for  flight  safety.  Developing  theoretical 
models  for  unsteady  flows  has  been  pursued  both  analytically  and  numerically.  But  because  of 
the  strong  nonlinear  character  of  the  governing  equations  (Navier-Stokes  equations),  one  still  must 
rely  on  experiments  to  understand  some  flow  characteristics  or  to  validate  theoretical  models.  The 
dramatic  advances  in  computer  power  in  the  past  few  years  have  made  Computational  Fluid  Dy- 
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namics  (CFD)  very  popular  among  engineers  and  researchers  to  overcome  some  of  the  difficulties 


of  solving  the  Navier-Stokes  equations. 

Flutter  is  a  phenomenon  generated  by  the  interaction  between  fluid  flow  and  structure. 
Therefore,  numerical  solutions  of  flutter  problems  require  the  use  of  CFD  coupled  with  Structural 
Dynamics  (SD)  solvers.  Analytical  solutions  of  the  unsteady  equations  are  available  only  for  a  few 
special  cases  with  little  practical  application.  Only  very  recently,  some  experimental  data  was  gen¬ 
erated  with  the  potential  for  validating  CFD/SD  codes,  for  example  the  data  presented  by  Schewe 
etal.  [31,44], 

The  study  of  flutter  started  because  of  the  earlier  development  of  fighter  aircraft,  especially 
during  World  War  I.  Some  designers  abandoned  the  classical  biplane  configuration,  with  its  inter¬ 
plane  bracing,  in  favor  of  the  mono-plane  wing  with  less  torsional  rigidity.  With  the  development  of 
faster  aircraft,  the  aeroelastic  behavior  of  the  airplane  became  more  important.  One  of  the  earliest 
books  published  on  this  subject  was  written  by  Bisplinghoff,  Ashley,  and  Halfman  [9]  in  1955.  It 
was  a  compilation  of  the  theoretical  and  experimental  tools  available  at  that  time  for  flutter  studies. 
The  analytical  models  were  mainly  based  on  the  small  disturbances  concept. 

The  first  significant  use  of  CFD  to  perform  flutter  calculations  is  due  to  Ballhaus  and  Goor- 
jian  [5]  in  the  late  1970’s.  They  studied  the  single-degree  pitch  stability  of  an  NACA  64A006  profile 
near  Moo  =  0.88.  The  main  concern  was  to  predict  correctly  the  motion  of  the  shock  with  oscillation 
of  the  profile,  particularly  the  phase  lag  with  respect  to  the  oscillatory  angle-of-attack. 

Since  then,  the  use  of  CFD/SD  solvers  to  analyze  flutter  problems  in  the  time  domain  be¬ 
came  more  and  more  frequent.  For  instance,  the  phenomenon  of  stall  flutter  was  studied  by  Eka- 
terinaris  and  Platzer  [17].  They  used  a  thin-layer  Navier-Stokes  solver  including  a  transition  model 
with  a  simplified  criterion  for  the  transition  onset.  The  numerical  solution  with  the  transition  model 
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showed  that  the  small  laminar/transitional  separation  bubble  forming  during  the  pitch-up  motion 


has  a  decisive  effect  on  the  near-wall  flow  and  the  development  of  unsteady  loads. 

A  time-domain  analysis  of  low-speed  airfoil  flutter  was  done  by  Jones  and  Platzer  [25]  in 
1996.  A  time-stepping  aeroelastic  code  was  used  to  perform  computations  for  two-airfoil  systems. 
The  configurations  studied  were  a  two-airfoil  tandem  and  an  airfoil  in  ground  effect.  Because  of 
the  excellent  agreement  with  past  results,  their  code  can  be  applied  as  a  feedback  loop  to  stabilize 
flutter  of  a  trailing  airfoil  actively,  to  investigate  wake  interference  in  rotary-wing  flowfields,  and  to 
simulate  flutter  in  ground  effect. 

Transonic  flutter  computations  were  the  subject  of  several  papers  published  by  Beran  et 
al.  [6,  10,  36].  The  main  interest  of  their  investigations  was  the  airfoil  flutter  boundaries.  They  used 
two  Euler  solvers  to  calculate  the  flutter  onset  for  a  pitch-and-plunge  airfoil  at  transonic  speeds. 
They  applied  the  Hopf-bifurcation  analysis  to  determine  the  flutter  onset  and  found  that  this  method 
is  precise  and  efficient  for  grids  typical  of  inviscid,  transonic  airfoil  calculations. 

Navier-Stokes  computations  for  flutter  of  an  airfoil  section  were  conducted  by  Weber  et 
al.  [54]  and  Weber  et  al.  [55].  In  the  former  work,  a  flutter  study  of  turbomachinery  blades  was 
conducted  by  using  a  full  Navier-Stokes  code.  In  the  latter,  Weber  et  al.  [55]  used  a  thin-layer 
Navier-Stokes  solver  to  perform  Limit  Cycle  Oscillation  (LCO)  computations  for  the  NLR  7301 
airfoil.  Comparisons  of  this  unbounded  flow  calculation  with  experimental  results  obtained  by 
means  of  a  wind  tunnel  test  suggested  significant  wall  interference  effects.  Castro  et  al.  [11]  found 
that  the  inclusion  of  the  wind  tunnel  walls  in  the  modeling  of  the  problem  can  significantly  improve 
the  agreement  of  numerical  and  experimental  results  for  steady  as  well  as  unsteady  computations. 

The  importance  of  tunnel  interference  is  well  recognized  [35],  but  reliable  quantitative  esti¬ 
mates  are  still  lacking.  For  instance,  Khalid  and  Mokry  [30]  have  presented  inviscid  calculations  for 
tunnel  wall  interference  effects  in  steady  flow  over  a  NACA  0012  airfoil  in  2-D  transonic  flow.  How- 
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ever,  no  Navier-Stokes  studies  seem  to  have  been  done  to  analyze  unsteady  transonic  interference 


effects. 


This  dissertation  is  a  continuation  of  the  work  presented  in  [11]  in  which  it  was  found  that 
both  the  porosity  and  the  way  of  applying  the  corresponding  boundary  condition  influenced  the 
results  of  both  steady  and  unsteady  computations  for  the  NLR  7301  airfoil  inside  a  wind  tunnel  at 
transonic  speeds.  In  that  work,  modeling  the  porosity  of  the  tunnel  wall  required  a  very  care ful 
construction  of  the  grid  in  which  the  cells  were  approximately  equally  spaced  at  the  wall  regions. 
This  was  necessary  because  some  grid  cells  were  treated  as  solid  walls  while  others  were  treated  as 
holes  in  order  to  achieve  the  desired  porosity.  For  example,  a  50%  porosity  was  modeled  by  treating 
four  consecutive  grid  cells  as  holes,  the  next  four  as  solid  walls,  and  so  on.  The  drawback  of  this 
approach  was  that  a  very  limited  number  of  possible  values  of  porosity  could  be  simulated  due  to 
the  finite  number  of  grid  cells  at  the  wind  tunnel  walls.  In  fact,  only  porosity  values  of  25%,  50%, 
and  75%  were  feasible.  Furthermore,  the  plenum  chamber  pressure  was  not  included  in  the  model. 

Comparisons  with  experimental  results  obtained  in  the  DLR-Gottingen  wind  tunnel  [31,  44] 
have  demonstrated  that  inclusion  of  wind-tunnel  effects  is  the  key  to  obtaining  better  agreement  of 
both  the  steady-state  pressure  distribution  and  the  flutter  characteristics  with  the  experiment. 

A  primary  goal  of  the  present  investigation  was  to  improve  the  porous  wind-tunnel  wall 
model  used  in  the  previous  study,  specifically,  providing  solution  uniqueness  and  including  the 
plenum  pressure.  As  earlier  mentioned,  the  previous  model  provided  a  number  of  ways  to  achieve 
the  same  porosity,  each  resulting  in  a  different  solution.  Furthermore,  it  was  desirable  to  elimi¬ 
nate  the  requirement  for  equally  spaced  gridding  along  the  tunnel  walls.  In  the  present  work,  the 
approach  presented  by  Mokry  et  al.  [35]  is  used.  The  normal  velocity  at  the  porous  wall  is  set 
proportional  to  the  pressure  difference  between  the  plenum  chamber  and  the  test  section  wall.  This 
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model  allows  a  much  more  continuous  variation  of  flow  parameters  and  does  not  require  an  equally 


spaced  grid  along  the  walls. 

In  this  dissertation,  the  effect  of  tunnel  wall  interference  on  transonic  flutter/limit-cycle 
prediction  is  investigated  in  further  detail.  To  this  end,  parametric  studies  of  porosity,  tunnel  height, 
and  Mach  number  are  conducted.  Numerical  solutions  are  obtained  for  an  airfoil  free  to  oscillate  in 
two-degrees-of-freedom  in  transonic  flow  and  the  results  are  compared  with  the  measurements  of 
Knipfer  et  al.  [31]. 

The  problem  of  transonic  flutter  of  an  airfoil  inside  a  wind  tunnel  is  solved  because  there  is 
neither  an  analytical  nor  an  experimental  method  of  eliminating  measurement  interference  caused 
by  wind  tunnel  walls  in  the  transonic  regime. 

B.  FLAPPING  AIRFOIL  STUDIES 

Flapping  wing  propulsion  has  been  used  by  birds  and  insects  for  many  millennia  but  barely 
by  mankind.  Although  the  concept  of  obtaining  thrust  by  means  of  a  plunging  airfoil  was  explained 
by  Knoller  [32]  in  1909  and  by  Betz  [7]  in  1912,  very  few  efforts  have  been  made  to  explore  or  to 
study  this  alternative  to  the  more  conventional  use  of  propellers. 

The  first  experimental  demonstration  of  the  feasibility  of  flapping  wing  propulsion  was  per¬ 
formed  by  Katzmayr  [29]  in  1922.  The  first  solution  for  incompressible  flow  past  flapping  airfoils 
was  presented  by  Birnbaum  [8]  in  1924.  Then,  in  1936,  Garrick  [19]  solved  the  problem  of  a 
sinusoidally  plunging  airfoil,  using  Theodorsen’s  theory  [48],  for  the  whole  range  of  reduced  fre¬ 
quencies.  This  solution  was  based  on  the  linear  theory  and  does  not  hold  for  detached  flows  typical 
of  high  reduced  frequencies. 
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One  of  the  earliest  numerical  models  for  a  3-D  flapping  wing  was  developed  by  Lan  [33] 


in  1979.  He  used  an  unsteady  quasi-vortex-lattice  to  calculate  the  propulsive  efficiency  and  thrust 
for  some  swept  and  rectangular  planforms  by  varying  the  phase  angles  between  the  pitching  and 
plunging  motions.  It  was  shown  that  aerodynamically  interacting  tandem  wings,  as  used  by  the 
dragonfly,  can  produce  high  thrust  with  high  efficiency  when  the  pitching  is  in  advance  of  the 
flapping  and  the  hind  wing  leads  the  fore  wing  with  some  optimum  phase  angle.  Linear  theory 
was  also  used  in  this  case,  therefore,  detached  flows  arc  not  captured  by  this  model. 

With  the  advent  of  CFD,  the  problem  with  detached  flows  could  be  overcome.  It  became 
possible  to  study  important  unsteady  phenomena,  such  as  the  dynamic  stall.  Among  the  recent  in¬ 
vestigations  done  on  this  subject,  the  work  of  Platzer  et  al.  [42,  43,  53]  deserves  special  attention. 
They  applied  a  thin-layer  Navier-Stokes  solver  improved  with  a  transition  model  based  on  Gostelow 
et  al.  [20]  and  determined  the  dynamic  stall  boundaries  for  the  NACA  0012  airfoil.  Furthermore, 
they  investigated  the  need  of  including  a  transition  model  into  the  CFD  solver  and  also  the  charac¬ 
teristics  of  dynamic  stall  and  flapping  airfoil  propulsion.  On  the  prediction  of  dynamic  stall,  Jones 
and  Platzer  [26]  showed  that  the  delay  in  the  dynamic  stall  onset  is  related  to  the  dynamic  pressure 
lag. 

Platzer  et  al.  [23,  28,  22,  51,  52]  also  studied  the  problem  of  flapping  wing  propulsion  in 
close  detail.  Numerically,  they  used  both  an  unsteady  panel  code  method  (UPOT)  coupled  with  a 
boundary-layer  solver  and  a  thin-layer  Navier-Stokes  solver  for  their  computations.  Experimentally, 
a  low-speed  wind  tunnel  and  also  a  water  tunnel  were  used.  Among  other  conclusions,  they  showed 
that  some  configurations  of  two  airfoils  can  produce  better  propulsive  efficiency  than  a  single  flap¬ 
ping  wing.  As  an  example  of  such  configurations,  the  flapping/stationary  airfoil  combination  in 
tandem  provided  a  40%  augmentation  in  propulsive  efficiency. 
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The  work  performed  by  Isogai  et  al.  [21]  is  also  worth  mentioning.  They  studied  the  effects 


of  dynamic  stall  on  propulsive  efficiency  and  thrust  of  a  flapping  NACA  0012.  It  was  determined 
that  a  phase  angle  between  the  pitch  and  plunge  motions  of  90  degrees  would  yield  optimum  values 
of  propulsive  efficiency  and  thrust.  A  similar  result  was  obtained  experimentally  by  Anderson  [3] 
and  Anderson  et  al.  [2].  Using  a  model  for  the  NACA  0012  immersed  in  a  water  tunnel,  it  was 
found  that  the  highest  propulsive  efficiencies  were  obtained  for  a  phase  angle  of  approximately  75 
degrees  between  pitch  and  plunge  and  motions  with  high  effective  angles  of  attack.  This  indicates 
that  the  presence  of  the  dynamic  stall  vortices  might  not  always  be  detrimental  to  thrust  generation. 

More  recently,  Dickinson  et  al.  [14]  experimentally  studied  the  characteristics  of  flight  of 
an  insect  wing.  It  was  found  that  the  most  important  mechanisms  in  insect  flight  are  what  they 
called  “rotational  circulation”  and  “wake  capture.”  In  the  latter  mechanism,  the  wing  benefits  from 
the  shed  vorticity  of  the  previous  stroke.  They  effectively  demonstrated  that  the  wing  can  produce 
positive  lift  by  interacting  with  its  own  wake. 

Ramamurti  and  Sandberg  [40]  also  studied  the  problem  of  an  oscillating  NACA  0012.  With 
a  finite  element  incompressible  flow  solver,  they  investigated  the  effects  of  varying  the  phase  angle 
between  pitch  and  plunge  motions.  Their  results  show  good  agreement  with  the  experiments  of 
Anderson  [3]  for  the  thrust  coefficient,  and  it  was  found  that  phase  angles  between  pitch  and  plunge 
around  90  degrees  deliver  the  best  propulsive  efficiency. 

Some  of  the  flapping  wing  configurations  investigated  in  recent  years  are  presented  in 
Fig.  1.1.  The  single  airfoil  (a)  has  been  the  subject  of  many  studies.  The  opposed-plunge  (c) 
combination  was  investigated  both  numerically  and  experimentally  by  Jones  and  Platzer  [28]  and 
Lund  et  al.  [34].  The  numerical  investigation  was  done  using  an  unsteady  potential  flow  code  mod¬ 
ified  to  compute  incompressible  flow  over  two  airfoils  oscillating  in  plunge  with  a  phase  angle  of 
180  degrees.  This  code,  called  USPOT,  was  originally  developed  by  Pang  [37]  in  1988.  It  is  im- 
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portant  to  recall  that  the  flow  solution  for  two  airfoils  arranged  in  a  biplane  fashion,  oscillating  in 


plunge  with  a  phase  angle  of  180  degrees,  is  equivalent  to  the  flow  over  a  single  airfoil  oscillating 
near  a  ground  plane.  This  remains  true  even  for  viscous  flow  solutions  as  long  as  an  inviscid  flow 
tangency  condition  is  enforced  along  the  ground  plane.  In  their  experimental  investigation,  Jones 
and  Platzer  [28]  built  and  tested  a  biplane  configuration.  Therefore,  the  single  airfoil  in  ground 
effect  solution  developed  in  this  dissertation  can  be  compared  with  their  experimental  results.  The 
experimental  paid  was  conducted  in  the  Naval  Postgraduate  School  low-speed  wind  tunnel. 

Apparently,  the  only  Navier-Stokes  computations  for  the  opposed-plunge,  or  biplane,  con¬ 
figuration  have  been  done  by  Tuncer  and  Kaya  [50].  Their  approach  consists  of  applying  C-grids 
for  each  airfoil  and  over-setting  them  with  a  Cartesian  grid  in  a  “Chimera”  fashion.  It  appeal's 
that  no  computations  were  done  for  such  a  configuration  using  a  deforming  grid  approach.  There¬ 
fore,  one  of  the  purposes  of  this  dissertation  is  to  analyze  this  problem  using  the  same  thin-layer 
Navier-Stokes  solver  used  for  investigating  the  transonic  flutter  of  an  airfoil  with  wind  tunnel  wall 
interference. 

The  problem  of  an  airfoil  oscillating  in  ground  effect  is  investigated  in  this  dissertation 
because  the  micro-air  vehicle  developed  at  the  Naval  Postgraduate  School  uses  a  biplane  wing  con¬ 
figuration  for  thrust  generation.  This  configuration  is  known  to  produce  more  thrust  per  wing  than 


the  single  wing  configuration  and  allows  a  much  better  mechanical  balance. 
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C.  COMMON  ISSUES 


The  flow  solver  and  the  aeroelastic  models  used  in  the  present  investigation  have  been 
tested  and  validated  extensively  in  previous  studies  for  a  variety  of  flow  conditions.  For  example, 
the  flow  solver  and  the  turbulence  models  have  been  tested  for  subsonic  flow  [43,  16,  18]  and  for 
transonic  flow  [15].  The  aeroelastic  model  has  been  implemented  and  tested  in  [27]  for  inviscid 
flow  calculations  and  in  [55]  for  viscous,  transonic  flow. 

The  transition  modeling  capability  of  the  code  was  demonstrated  in  [55,  43].  Nevertheless, 
because  previous  results  of  [55]  showed  that  transition  only  slightly  improved  the  numerical  results, 
transition  modeling  was  not  used  in  the  present  investigation,  assuming  fully  turbulent  flow. 

The  influence  of  the  different  turbulence  models  available  for  the  present  code  in  the  numer¬ 
ical  predictions  of  transonic  flutter  was  also  studied  by  Weber  et  al.  [55].  It  is  known  that  turbulence 
models  are  calibrated  only  for  steady-state  flows.  Nevertheless,  it  is  common  practice  to  apply 
them  for  unsteady  flows  without  any  further  modification.  It  was  found  during  this  work  that,  for 
one-equation  turbulence  models  applied  to  moving  grid  problems,  some  modifications  in  their  im¬ 
plementation  is  necessary.  Specifically,  the  time  metrics  should  be  included  in  the  computation  of 
the  contravariant  velocity  used  in  the  convective  terms,  and  a  time-accurate  integration  of  the  trans¬ 
port  equation  should  be  performed.  This  modification  was  discussed  with  Dr.  Steven  Allmaras  [1] 
and  represents  one  of  the  contributions  of  the  present  investigation. 

The  ability  of  the  multi-block  version  of  the  flow  solver  to  predict  accurately  the  flow  over  a 
stationary  airfoil  in  a  wind  tunnel  by  including  porous  wall  effects  was  demonstrated  in  the  previous 
work  [1 1]. 

The  problems  of  flutter  of  an  airfoil  inside  a  wind  tunnel  and  of  an  oscillating  airfoil  near 
a  wall  can  be  numerically  solved  in  a  very  similar  manner.  For  the  former,  two  wall  boundary 
conditions  should  be  applied  both  above  and  below  the  airfoil.  For  the  latter,  only  one  wall  boundary 
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Figure  1.2.  Comparison  of  the  Configurations 


condition  is  necessary  below  the  airfoil.  For  the  upper  boundary,  free-stream  boundary  conditions 
arc  applied.  Therefore,  the  same  solver,  only  with  minor  changes  owing  to  specific  characteristics 
of  each  configuration,  can  be  used  for  numerical  computations  of  both  problems. 

Time-accurate  unsteady  computations  require  large  computing  times  until  a  reasonable 
number  of  cycles  is  reached.  Computational  efficiency  can  be  improved  by  taking  advantage  of 
parallel  architectures.  The  computational  domain  is  divided  into  three  blocks  to  enable  discretiza¬ 
tion  of  the  full  geometry,  for  both  problems  under  investigation.  This  makes  the  approach  suitable 
for  parallel  computation.  Therefore,  a  newly  developed  parallel  version  of  the  code  uses  three 
Pentium  II  400  MHz  PC’s  to  carry  out  the  computations,  assigning  one  block  to  each  processor. 
The  boundary  conditions  arc  transferred  from  one  machine  to  another  via  a  specific  library  called 
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Message  Passing  Interface  (MPI).  The  cluster  of  Linux  PC’s  is  effective  in  performing  the  parallel 


computations  and  is  also  effective  in  reducing  the  wall-clock  time. 


D.  DISSERTATION  OVERVIEW 

The  theoretical  background  for  the  computational  codes  used  in  this  work  is  presented  in 
Chapter  II.  The  thin-layer  Navier-Stokes  equations  in  conservative,  unsteady,  generalized  coordi¬ 
nate,  and  nondimensional  form  arc  discussed  along  with  some  important  issues,  such  as  the  numeri¬ 
cal  scheme,  turbulence  models,  deforming  grids,  and  transition  models.  The  necessary  adjustments 
for  working  with  moving  and  changing  grids  arc  pointed  out.  A  brief  explanation  of  the  two-degree- 
of-freedom  structural  dynamics  model  is  also  given  in  Chapter  II. 

Some  more  specific  topics  concerning  the  present  work  are  presented  in  Chapter  III.  It  con¬ 
tains  an  explanation  of  why  a  three-block  grid  configuration  is  used  and  how  a  parallelized  version 
of  the  program  is  implemented.  The  treatment  of  boundary  conditions,  one  of  the  most  impor¬ 
tant  issues  concerning  unsteady  aerodynamics,  is  also  discussed  in  Chapter  III.  Specific  boundary 
conditions  for  both  the  transonic  flutter  in  the  wind  tunnel  and  the  airfoil  in  ground  effect  configu¬ 
ration  are  explained  in  more  detail.  Furthermore,  a  combined  pitch/plunge  motion  is  imposed  to  a 
NACA  0012  airfoil  in  order  to  evaluate  the  performance  of  two  different  boundary  conditions  at  the 
trailing  edge. 

The  validations  of  the  computational  code  used  in  this  work  arc  presented  in  Chapter  IV, 
specifically  with  respect  to  the  moving  or  deforming  grid  capability  of  the  solver.  A  transonic 
steady-state  computation  is  performed  to  show  that  the  program  does  not  generate  large  deviations 
from  the  solution  when  moving  grid  points  arc  crossing  shock  waves.  A  subsonic  numerical  test  is 
conducted  to  show  that  a  multi-grid  discretization  yields  a  similar  solution  as  a  single-grid  approach. 
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At  the  end  of  this  chapter,  a  comparison  with  a  potential  flow  solver  solution  is  done  for  a  biplane. 


opposed-plunge  configuration.  The  purpose  of  this  chapter  is  to  show  that  the  theory  presented  in 
Chapters  II  and  III  is  adequate  for  treating  deforming  grid  problems. 

The  results  obtained  for  the  two  problems  under  investigation  arc  presented  in  Chapter  V. 
In  the  first  section  of  this  chapter,  the  details  about  the  results  for  the  flapping  airfoil  in  ground 
effect  configuration  arc  discussed.  A  NACA  0014  airfoil  is  allowed  to  oscillate  in  the  pure  plunge 
mode  near  a  wall.  The  influence  of  some  parameters,  such  as  the  Reynolds  and  Mach  numbers  is 
investigated.  Fully  turbulent  computations  are  done  for  the  Spalart-Allmaras  and  Baldwin-Lomax 
turbulence  models.  A  fully  laminar  simulation  is  performed  for  a  Reynolds  equal  to  10,000  as  well. 
The  second  part  of  this  chapter  is  dedicated  to  the  transonic  flutter  with  tunnel  wall  interference 
of  the  NLR  7301  airfoil.  The  influence  of  the  porous  wall  boundary  condition  is  studied  in  close 
detail.  The  porosity  parameter  is  found  to  change  the  steady-state  solutions  significantly.  Therefore, 
porosity  parameters  influence  on  the  limit  cycle  oscillations  is  also  evaluated.  Because  the  Mach 
number  is  close  to  the  transonic  dip,  this  parameter  is  varied  in  order  to  evaluate  the  changes  in 
the  flutter  characteristics  of  the  NLR  7301  section.  The  height  of  the  wind  tunnel  test  section  is 
increased  to  assess  correlations  of  wind  tunnel  tests  with  unbounded  flow  situations. 

Finally,  the  conclusions  of  the  present  work  for  the  two  problems  under  investigation  and 
the  recommendations  for  future  studies  arc  discussed  in  Chapter  VI. 
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II.  THEORETICAL  BACKGROUND 


A.  UNSTEADY  AERODYNAMICS 

1.  Governing  Equations 

The  solution  for  the  problem  of  an  unsteady,  compressible,  and  viscous  flow  of  a  Newtonian 
fluid  is  obtained  by  solving  the  Navier-Stokes  (N-S)  equations.  They  are  presented  in  matrix,  non- 
dimensional,  and  Cartesian  coordinate  form  in  Eq.  (2.1): 


dQ  dF  dG  ,  fdR  dS\ 

!k+!h+!k~Re  + 


(2.1) 


where  Q  =  (p,pH,pw,e)r  is  the  dependent  variable  and  represents  a  vector  which  components  are 
the  flow  state  variables. 

The  reference  values  for  non-dimensionalization  are  the  chord  length  c,  the  free  stream 
density  poo,  the  free  stream  speed-of-sound  a„,  the  time  c/a„,  and  the  specific  energy  p^rQ. 

The  fluxes  F,  G,  R ,  and  S  are  given  by: 
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'  'N 

'  \ 

P  u 

p  w 

p  u2  +  p 

>3  G  =  < 

p  uw 

p  uw 

p  W2  -1-  p 

(e  +  p)u 

(e  +  p)w 

/ 

(2.2) 


0 

0 

R=  < 

^ XX 

>  ,  S  =  < 

Ur 

^xz 

x zz 

^XxU  +  T'ZX^  Qx 

V.  / 

TXZU  +  ^zzw  ~  Qz 

K 

where: 


X-«  ~^(tx  + 


(2.3) 


q.x  =  -fc 


"52  5 


qz=-k^,  X  —  —2/7/3  . 


Pressure  is  related  to  the  other  variables  through  the  equation  of  state  for  an  ideal  gas 


p  =  p RT  . 


(2.4) 


The  N-S  equations  have  analytical  solutions  only  for  a  few  simple  problems.  Usually, 
their  solution  involves  a  numerical  procedure  using  a  digital  computer.  This  technique  is  called 
Computational  Fluid  Dynamics  (CFD)  and  can  be  performed  either  with  the  Finite-volume  method 
or  the  Finite-difference  method.  Both  methods  require  the  generation  of  a  computational  grid  to 
discretize  the  flow-field.  The  grid  can  be  adaptive  to  the  boundaries  of  the  domain  where  the  N-S 
equations  arc  being  solved.  Such  a  situation  demands  the  N-S  equations  to  be  written  in  generalized 
coordinates  (x,  q,  and  Q: 
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Qx  +  Fk  +  G^  =  Re-l{R%  +  \)  , 


(2.5) 


where: 


Q  =  J  lQ 

j  =  (x^z^-x^y1 

F  =  i,Q  +  ixF  +  izG 

R  =  lxR  +  lzS 

g  =  Iq+\xf  +  Izg 

S  =  t:R  +  & 

It  =J~%=x^Zt-xxz^ 

lt  =  J-l£>t  =  xzZi~-xt:Zz 

4v  =  J~% X  =  ^ 

L  =  j-%  = 

4z  =  7  %  =  -X$ 

Zx  =  J-%  =  x$ 

(2.6) 


(2.7) 


Equation  (2.5)  represents  the  Navier-Stokes  equations  in  conservative,  non-dimensional 
form  and  is  written  in  terms  of  the  computational  domain  variables.  They  are  also  valid  for  moving 
or  deforming  grids  because  the  Jacobian,  J,  is  considered  to  be  a  function  of  time  (see  Appendix  A 
for  details). 

A  common  simplification  to  the  Navier-Stokes  equations  is  the  thin-layer  approximation 
where  the  term  R%  on  the  RHS  of  Eq.  (2.5)  is  neglected  in  comparison  with  Sr- 


dzQ  +  d^F  +  d^G  =  Re~ld^S  . 


(2.8) 


Substituting  Eq.  (2.2)  into  Eq.  (2.6): 


17 


Q  =  J  \ 


P 

pu 
p  w 


\ 

P  u 

>  F  =  J-l< 

puU  +  qxp 

pwU  +  t,zp 

{e  +  p)U-^tp 

/ 

f  N 

pVT 

0 

A  1 

G  =  J-[  < 

puW  +  C,xp 

A  1 

>  s  =  y-J< 

pmiu^  +  (p/3)m2C,x 

pwW  +  t,zp 

pmiw^  +  {p/3)m2C,z 

(e  +  p)W  -t^p 

V  / 

pmitni  +  (p/3)m2m4 

\  / 

where: 


(2.9) 


m  i  =  {,1+Q 
m2  =  'QxU'q  +  t,zw^ 

m3  =  3^(m2  +  w2)/2+  [(y—  l)Pr]_13^(fl2) 
m4  =  C,xu  +  £zw 

The  terms  t/  and  IT  are  the  contravariant  velocity  components  given  by: 

U  =  +  ^A-w  +  £,zw 

W  =  C,r  +  C; rH  + 

Pressure  is  related  to  the  other  variables  through  the  equation  of  state  for  an  ideal  gas: 


(2.10) 


(2.11) 


P={ Y-  1)  [e-p{u2  +  w2)/2] 


(2.12) 


2.  Turbulence  Models 

Turbulence  models  arc  introduced  to  solve  the  closure  problem  generated  by  time  averaging 
the  N-S  equations.  They  arc  divided  into  categories,  such  as  Algebraic,  One-Equation,  and  Two- 
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Equation  models.  The  method  used  in  this  work  is  implemented  with  one  Algebraic  turbulence 
model  called  Baldwin-Lomax  (B-L)  and  one  One-Equation  model  called  Spalart-Allmaras  (S-A). 

The  effects  of  turbulence  arc  translated  in  terms  of  an  eddy  viscosity  coefficient  (pf)  and 
new  heat  flux  terms.  The  fluid  viscosity  p  is  replaced  by  p/  +  p, ,  where  p/  represents  the  laminar 
viscosity,  and  k/cp  =  p/ Pr  by  k/cp  =  p/  / Pr  +  p,  / Pr, . 

a.  Baldwin-Lomax  Turbulence  Model 

The  Baldwin-Lomax  turbulence  model,  presented  in  [4],  represents  a  two-layer 
algebraic  eddy  viscosity  model  in  which  the  term  p,  is  described  by: 


= 


(P/ )  inner-,  )  y crossover 

{fit) outer )  y  A  y crossover 


(2.13) 


where  y  is  the  normal  distance  from  the  wall  and  y crossover  is  the  smallest  value  of  y 
at  which  values  from  the  inner  and  outer  equations  are  equal. 

The  equations  for  the  inner  sub-layer  are  based  on  a  Prandtl-Van  Driest  formulation: 


{fit)  inner  —  pA"|to| 


where: 


l  =  Ky[l  —  e  y+  /A+] 

M  =  \/(£-i)2 

_  pwuxy VP w^wy 

y  ~  /jw  ~  /iw 

The  following  equations  are  used  for  the  outer  region: 


{fit)  outer  —  KCcPpFwakeFkleb{y) 


(2.14) 


(2.15) 


(2.16) 
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where  K  is  the  Clauser  constant,  Ccp  is  an  additional  constant,  and 


F wake  —  knin(ymaxF nax  iCwkymaxUrfif/  F wax) 

F(y)  =  yM[i-e_yf/A+] 


1  +  5.5 

y  jmax  J 


n  -i 


(2.17) 


Fklebiy) 

Udif  =  (y/u2  +  w2)  max  (V  u2  +  w2),. 


The  variables  ymax  and  Fmax  are  determined  from  the  function  F  (y) .  In  wakes, 
the  exponential  term  of  this  function  is  set  equal  to  zero.  The  function  Fkieb{y)  is  the  Klebanoff 
intermittency  factor.  Udif  is  the  difference  between  maximum  and  minimum  total  velocities  in  the 
profile  for  a  fixed  x  station.  The  second  (min)  term  in  the  equation  for  Udif  is  taken  as  zero  for 
wakes.  These  equations  avoid  the  necessity  of  finding  the  outer  edge  of  the  boundary  layer. 

The  constants  in  the  equations  for  the  B-L  model  arc  given  below: 


A+=  26  K  =  0.4 

Ccp  =1.6  K  =  0.0168 

Ckieb  =  0-3  Pr  =  0.12 

Cwk  =  0.25  Pr >  =  0.9 

b.  Spalart-Allmaras  Turbulence  Model 


(2.18) 


The  S-A  turbulence  model  is  suggested  in  [45].  Reynolds-averaged  Navier-Stokes 
equations  and  a  transport  equation  arc  solved  for  this  turbulence  model.  The  formulation  for  the 
eddy  viscosity  v,  is  given  by: 

(2.19) 

vf  =  v/vi  fv  t  =  1  =  1 

l  TSl 


UjUj  —  2vtSi  j  Sij  —  ^ 


dUi  ,  Mj 
dxj  '  dxi 
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where  v  is  the  molecular  viscosity,  v  is  the  working  variable  and  obeys  the  transport  equation  for 


this  model: 


§  =  CM[l-/r2]5v  +  i[V-((v  +  v)Vv)+C,2(Vv)2] 

-[cwxfw-c-$fa]  [^]2  +  fnAU2 


(2.20) 


and 


^s+^fv2  /v2  =  i-t^ 


fw  —  g 


l  +  c' 


U6+4? 


.6  1  t/6  , 

—  1  g=  r  +  cw2{t  -  r)  r  = 


(2.21) 


§K2d2 


S  represents  the  magnitude  of  the  vorticity  and  cl  is  the  distance  to  the  closest  wall.  The  boundary 
condition  for  the  wall  is  obtained  by  setting  v  =  0.  The  functions  fa  and  f,  \  are  expressed  by: 


ft 2  =  ct3  exp{-ctAl2 ) 

fn  =  Qi  gt  exp  (-c,2-$p\cl2  +  g2d2]}  (2-22) 

gt  =  min(0.l,AU  /w,Ax) 

where  dt  is  the  distance  from  the  field  point  to  the  trip,  which  is  on  the  wall,  0),  is  the  wall  vorticity 
at  the  trip,  AU  is  the  difference  between  the  velocity  at  the  field  point  and  the  velocity  at  the  trip, 
and  Ax  is  the  grid  spacing  along  the  wall  at  the  trip. 

The  constants  used  for  this  model  arc 


q,i  =  0.1355 

cw\  —  Cbl/K2  +  (1  +  Cb2)/G 

C/1 

=  1.0 

Cb2  =  0.622 

Cw  2  =  0.3 

C/2 

=  2.0 

a  =  2/3 

cw3  =  2.0 

C/3 

=  1.1 

K  =  0.41 

Cvi  =  7.1 

C/4 

=  2.0 

(2.23) 


The  constants  used  in  this  turbulence  model  are  calibrated  for  steady  flows.  In 
this  work,  the  constants  arc  assumed  to  be  good  for  unsteady  flows  as  well.  This  is  a  necessary 
assumption  because  there  arc  no  calibrations  for  unsteady  flows  so  far.  More  importantly,  the  coding 
of  this  turbulence  model  assumes  that  the  grid  is  not  changing,  i.e.,  the  grid  is  not  a  function  of  time. 
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In  other  words,  this  turbulence  model,  as  it  is  usually  coded,  is  calibrated  for  steady  flows  and  is 


not  changing-grid  compatible.  The  latter  is  because  the  contravariant  velocities  arc  coded,  in  the 
turbulence  model  subroutine,  simply  as: 

U  =  t,xu  +  t,zw  and  W  =  L,xu  +  (2.24) 

For  the  turbulence  model  to  be  moving-  or  deforming-grid  compatible,  it  must  use  the  following 
equation  for  the  contravariant  velocities: 

U  =  %t  +  ^xu  +  t,zw  and  W  =  C,t  +  C,xu  +  t,zw  (2.25) 

The  terms  ^  and  C,t  represent  changes  in  the  grid-conforming  coordinates  with  time  and  arc  not 
related  to  the  fact  that  the  flow  is  unsteady.  There  might  be  a  situation  in  which  the  flow  is  steady 
and  the  grid  is  changing,  such  as  in  an  adaptive  grid  solution,  and  those  time  metrics  should  be 
included  in  the  equations  of  the  turbulence  models. 

In  a  personal  communication  with  Dr.  Steven  Allmaras  [1],  he  confirmed  that  the 
changes  made  to  the  conservation  equations  (mass,  momentum,  energy),  to  incorporate  unsteady 
and  moving  grids,  should  be  duplicated  for  the  S-A  model.  Another  issue  mentioned  by  Dr.  All¬ 
maras  was  the  time  integration.  Most  turbulent  codes,  as  well  as  the  present  one,  use  a  decoupled 
relaxation  (or  time  integration)  between  the  conservation  equations  and  the  turbulence  model.  This 
decoupled  time  integration  introduces  a  first-order  error  (O(At))  into  the  results.  This  error  can 
become  significant  for  high  frequencies.  A  solution  for  this  problem  would  be  implementing  a 
subiteration  process  at  each  time  step. 

3.  Transition  Models 

The  transition  from  laminar  to  turbulent  flow  is  a  highly  important  issue  on  flow  calculations 
because  it  has  a  strong  influence  on  the  predictions  of  losses,  separation,  and  stall.  Transition  also 
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has  a  significant  effect  on  the  computed  values  of  wall  shear  stresses  and  surface  heat  transfer. 
Transition  criteria  used  today  in  CFD  codes  are  based  mainly  on  the  work  of  Gostelow  et  al.  [20]. 

Transition  models  start  with  the  prediction  of  the  transition  onset  location  and,  then,  move 
on  to  the  evaluation  of  the  transition  length.  In  this  work,  transition  onset  location  is  predicted 
using  Michel’s  Criterion  [12]  and  the  transition  length  is  evaluated  based  on  the  model  proposed  by 
Gostelow  et  al.  [20]. 

a.  Michel’s  Criterion 

This  criterion  assumes  that  transition  is  initiated  when  the  Reynolds  number  based 
on  the  momentum  thickness.  Re q,  and  the  Reynolds  number  based  on  x,  Rex,  satisfy  the  equation: 

Ree  =  1 . 174  (l  +  22’400  |  Re  °-46  where  ReQ  =  UeQ /v„  .  (2.26) 

V  ReXtr  J 

Equation  (2.26)  is  solved  for  xtr  and  this  value  is  used  as  x,  in  the  Gostelow  et 

al.  [20]  model. 

b.  Gostelow  et  al. 

The  idea  behind  this  model  is  that  the  turbulent  viscosity  should  be  weighted  by  an 
intermittency  function  to  provide  a  smooth  transition  from  laminar  to  turbulent  flow.  This  is  done 
according  to  Eq.  (2.27): 

Rtrans  —  Y (x^Rturb  •  (2.27) 

The  intermittency  function  in  the  transitional  region  is  given  by 

rx  (j  f dx\  fx 

y(x)  =  1  -exp  )  l  tanedxj  (2.28) 

where  a  is  the  spot  propagation  parameter,  e  is  the  spot  spreading  half  angle,  and  A.e  is  the  pressure 
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10.86  x  1(T3  exp (2.134Xe,  ln(^)  -  59.23X0,  -0.564 Info)),  for)i0,  <0 

(2.33) 

N(le  =  0)  x  exp(-10v/X^’),  for  A,0f  >  0 

and  where  q,  denotes  the  free-stream  turbulence. 

The  spot-propagation-rate  and  the  spot  spreading  half-angle  asymptotically  ap¬ 
proach  a  maximum  value  for  high  negative  values  of  A.0,  but  n  is  allowed  to  increase  to  infinity 
for  high  negative  values  of  X0f ,  where  L0,  is  the  pressure  gradient  at  the  transition  onset  location,  xt. 

The  value  of  the  intermittency  parameter,  y(x),  is  zero  for  x  <  xt,  and  increases 
downstream  from  the  transition  point,  asymptotically  to  a  maximum  value  of  unity,  which  corre¬ 
sponds  to  fully-turbulent  flow.  The  value  of  x,  is  usually  given  by  Michel’s  criterion. 


4.  Numerical  Technique 

The  Thin-Layer  Navier-Stokes  (TLNS)  equations  are  discretized  using  an  alternate  direc¬ 
tion  implicit  (ADI),  third-order  accurate  in  space,  second-order  accurate  in  time,  finite-volume 
scheme,  which  can  be  represented  by: 
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(2.34) 


/  +  /.4(V^++A?ir1)/>]AQ*,  =  -^ 

'+MvA!i ■,»)<■]  (eft1  -ey  =  aq*, 

where: 

^  =  (e^-ey +^+1/2,*- Aw)' 

+^(<^,£+1/2  -  G,',A-1/2)P  -  ^-1/zC(%+1/2  -  ^i,k— l/l)^ 

A  A  A  A  A  /V 

The  variables  A,  B,  and  M  are  the  flux  Jacobian  matrices  and  are  defined  as  A  =  0F/0Q, 
B  =  dG/dQ ,  and  M  =  dS/dQ,  respectively.  A  flux  splitting,  according  to  Steger  and  Warming  [46], 

A  A  A  A  |  A  A  A  j  A 

is  applied  to  matrices  A  and  B,  where  A  =  A+  +  ,4  and  B  =  B  \-  B  . 

The  h  quantities  are  defined  as  hr  =  Ax /Ac,  and  hr  =  Ax/A^.  V,  A,  and  5  are  the  forward, 
backward,  and  central  difference  operators,  respectively. 

The  variables  Fi+  \/2^,  Gt^+ 1/2,  and  S^k+ 1/2  are  numerical  fluxes.  The  superscript  (.)"  de¬ 
notes  the  physical  time  step  and  the  superscript  {.)p  is  related  to  Newton  sub-iterations  within  each 
physical  time  step,  which  are  used  to  improve  time  accuracy.  These  sub-iterations  minimize  the 
linearization  and  factorization  errors  and  help  drive  the  left-hand  side  of  Eq.  (2.34)  to  zero. 

Inviscid  numerical  fluxes,  Fi+i/2^  and  G,-^+ 1/2,  are  evaluated  by  means  of  the  Osher’s  third- 
order  accurate,  upwind-biased  scheme  [13].  Linearization  of  the  left-hand  side  of  Eq.  (2.34)  is 
performed  by  evaluating  the  flux  Jacobian  matrices,  A  and  B,  with  the  Steger  and  Warming  flux- 
vector  splitting  [46].  The  viscous  numerical  flux  Sj^+  \/2  is  computed  with  second-order  central 
differences.  Furthermore,  a  Total  Variation  Diminishing  (TVD)  flux  limiter  suggested  by  Rai  and 
Chakravarthy  [39]  is  applied  to  minimize  numerical  oscillations  at  shocks  developed  at  transonic 
speeds.  Grid  movement  and  deformation  compatibility  is  obtained  by  using  the  methodology  pre¬ 
sented  by  Tamura  and  Fugii  [47], 
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a.  Steger  and  Warming  Flux  Vector  Splitting 

The  flux  vectors  F(Q)  and  G(Q )  can  be  split  into  two  parts.  This  is  accomplished 
by  considering  the  fact  that  the  Jacobian  matrices  associated  with  the  fluxes  have  a  complete  set  of 
lineal-  independent  eigenvectors  and,  then,  can  be  split  into  subvectors.  The  equations  for  the  flux  F 
are  given  below  and  the  ones  for  flux  G  are  similar. 


where: 


where: 


A,  C7  JT  A  A  A  1  A, 

F=-^Q  =  AQ  =  LAL  lQ: 
dQ 


0 

0 

0 

u 

0 

0 

0 

0 

X2 

0 

0 

0 

u 

0 

0 

0 

0 

X3 

0 

0 

0 

u  +  c 

0 

0 

0 

0 

X4 

0 

0 

0 

u  —  c 

Any  eigenvalue  A,/  can  be  written  as: 

X,  =  Xf  +  Xj 


1+  _  ^  +  N 

Ai  ~  o 


,  and  V  = 

2  '  2 


Therefore,  matrix  A  can  be  split  into  to  other  diagonal  matrices, 


A  =  A+  +  A~ 


(2.36) 


(2.37) 


(2.38) 


(2.39) 


(2.40) 


where  the  diagonal  entries  of  A+  and  A  arc  and  ,  respectively. 

Finally,  applying  Eq.  (2.40)  to  Eq.  (2.36),  the  numerical  flux  F  can  be  expressed 
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as: 


F  =  £(A++A-)£-1(2 
F=(A++A-)<2 


(2.41) 


where: 


A+  =  IA+I^1,  A~  =  £A~£_1  ,  and  A=A+AAT  (2.42) 


b.  Osher  Upwind  Scheme 


A  detailed  explanation  of  the  Osher  scheme  is  presented  in  [13].  Nevertheless, 
some  comments  must  be  made  to  better  explain  the  interaction  of  the  Osher  method  with  the 
methodology  for  treating  deforming  grids.  In  order  to  facilitate  notation,  the  hats  on  the  fluxes 
will  be  dropped. 

For  a  first-order  Osher  scheme,  the  numerical  fluxes  are  updated  using: 


<2',+1  =  £;•'■ 


Ax 

A^ 


rQi  ,  rQi+i 

/  A+dQ+  /  A~dQ 
J  Qi-  i  0 


(2.43) 


ing  to: 


The  integrals  of  Eq.  (2.43)  are  evaluated  on  the  subpaths  shown  in  Fig.  2.1  accord¬ 


s£ _A+dQ  +  f°;+lA  dQ=  F{^Qi+x)-F{^Qi) 

+ jS  ,-'•  (f )  1  ds+ !%;:%. A+  (f  )2ds 
+  ^a-i/3A+  (d^)3ds~fQi+  7  A+  {^)4ds 

{^f)5ds-&,3A+  (* )6ds 


(2.44) 


The  point  where  the  metrics  are  evaluated  can  be  chosen.  For  instance,  one  can 


define: 
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i-2/3 


i- 1/3 


i+1/3 


i+2/3 


i-1  i  i+1 


Figure  2.1.  Paths  of  Integration  for  the  Osher  Scheme 


F{£>,  2 )  =  F(ki  1/2;  2 )  when  evaluating  /|'  A+dQ 

(2.45) 

^(£,2)  =^(^+1/2,2)  when  evaluating  /J'+1^~^2 

The  choice  of  the  mid-point  of  the  interval,  ^i±\/2,  is  quite  convenient  because  it 
requires  only  about  half  of  the  computations  for  each  time  step.  This  is  possible  because  some 
terms  that  are  evaluated  in  the  calculation  of  the  numerical  flux,  for  a  particular  interval,  arc  equal 
to  terms  evaluated  for  the  previous  one.  Thus,  these  terms  do  not  need  to  be  evaluated  again, 
saving  computational  time.  This  is  called  telescopic  property  and  is  explained  in  more  detail  in 
[13].  Another  consequence  of  using  mid-point  metrics  is  that  it  also  allows  the  numerical  scheme  to 
satisfy  the  conservation  in  the  discretized  form  of  the  equations,  which  is  one  of  the  requirements 
for  a  moving  grid  solution. 

For  a  third-order  accurate  Osher  scheme,  Eq.  (2.43)  is  changed  to: 


ef-e"  ,  f+i/2  1/2 


-  4^  [AF+  (G,-_1 ,  Qu  5i+1/2)  -  AF+  (0;_2,  2.  1, 5,-—!/ 2)]p 
+  ^[AF-(Qi+hQi+2^i+li2)-AF-(QhQi+h^-l/2)}p 


(2.46) 


TVD  Flux  Limiter 


The  TVD  scheme  used  in  the  present  work  is  based  on  [39].  The  fluxes  AF+  and 
A F  on  Eq.  (2.46)  are  modified  to  achieve  the  TVD  property.  The  flux  limiters  arc  obtained  with 
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the  help  of  an  entropy  function,  defined  as: 


V(Q)  =  -plogip/p*)  (2.47) 

The  idea  behind  this  flux  limiter  is  to  change  the  fluxes  when  the  gradient  of  this 
entropy  function  exceeds  a  certain  limit,  which  is  what  occurs  when  a  shock  wave  exists.  Then,  the 
gradient  of  the  entropy  function  described  by  Eq.  (2.47)  can  be  expressed  as: 

VQ  =  S/QV  =  -^  p1^  e+^-^-Y-l  +  (og^,-pn,-pv,p  (2.48) 

The  modification  of  the  fluxes  is  done  for  each  subpath.  Defining  DVq"'1  as  the 
change  in  the  entropy  function  along  he  subpath  m, 

DV^  =  (Ve)i+IB/3  -  (VQ)i+{m- 1)/3  ,  m  =  1,2,3  (2.49) 

and  the  modified  fluxes  (AF+)(m)  and  (AF~)(m)  are  given  by: 

[AF+(G/-1,G/,y](m)  =  [AF+(G/,Gi+1,y](M)  X  max[0,min(^N/D)\  (2.50) 

where 

1-0  <  (|>  <  2.0 

N=<  [AF+(Q/_1,Q,^)]W,[Dye(G;,Gi+1,^)]W  >  (2-51) 

D=<  [AF+(QhQi+iM{m\[DVQ(QhQi+  i,5)](m)  > 
and 

[AF-(G«,  G/+  i,5)]W  =  [AF-(Gi^  G/,5)]W  x  max[0,min(^N/D)\  (2.52) 

where 

1-0  <  <>  <  2.0 

N=<  [A F-{QuQi+lM{m\[DVQ{Qi-uQiM{m)  >  (2‘53) 

D=<  [A F-{Qi^Qi^m\[DVQ{Qi^iM{m)  > 

The  variables  N  and  D  represent  the  inner  product  of  the  vectors  contained  within  the  symbols  <  >. 
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d.  Deforming-Grid  Compatibility 

The  N-S  equations  represented  by  Eq.  (2.8)  are  compatible  with  deforming-  and 
moving-grid  problems,  as  shown  in  Appendix  A.  Theoretically,  the  N-S  equations  can  handle  prob¬ 
lems  in  which  both  q,  and  'Q,  are  different  from  zero.  Computationally,  the  numerical  scheme  must 
be  conservative  both  in  time  and  space  in  order  to  guarantee  deforming-grid  compatibility. 

Thomas  and  Lombard  [49]  suggest  the  use  of  the  Geometric  Conservation  Law 
(GCL)  for  solving  problems  involving  a  changing  or  deforming  grid.  As  shown  in  Appendix  A,  the 
GCL  is  analytically  satisfied  by  the  unsteady  N-S  equations  written  in  generalized  coordinates  form. 
Nonetheless,  in  order  to  satisfy  the  GCL  numerically,  the  metrics  associated  with  computations  of 
the  numerical  fluxes  must  be  calculated  using  a  conservative  scheme.  They  suggest  the  following 
conservative  scheme  for  the  3-D  metrics: 


lx  =  Ovk  -  Wit  =  Wk  ~  M? 

(2.54) 

lx  =  (^)n  -  Ovh  etc- 

Lor  2-D  metrics,  Eq.  (2.54)  is  reduced  to  the  form  presented  in  Eq.  (2.7)  and  the 
scheme  is  automatically  conservative  in  space.  Therefore,  the  numerical  scheme  used  in  the  present 
work  is  conservative  in  space.  Thomas  and  Lombard  [49]  also  suggest  the  inclusion  of  an  extra  term 
in  the  discretized  equations  to  achieve  moving-grid  compatibility.  According  to  their  approach,  the 
Osher  scheme  presented  in  section  b  would  become: 


U- '  +  +  a5a-,)"]  A e*,  = 

+  (Q%1  -  Q^)  =  IJ-'AQ 


Lk 


(2.55) 
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where: 


Now, 

A{J-1 


=  J-HQU  ~  QW  +  MJ~X)Qb  +  h(Fi+l/2,k  ~  Pi- 1/2*)' 

(2.56) 

+^{Gi,k+\/2  ~  Gi,k-i/2)p  ~  Re^h^{SiMi/2  ~  Si,k- itiY 

the  (7)  terms  represent  Jacobians  with  respect  to  Q,  or  d(.)/dQ,  instead  of  Q  and  the  term 
)  is  obtained  by  the  GCL  (see  Appendix  B): 


A(/-1)  =  -At[(|^  +  (^] 


(2.57) 


One  of  the  differences  between  Eq.  (2.55)  and  Eq.  (2.34)  is  the  way  they  arc  solved. 
The  former  solves  for  Q  while  the  latter  solves  for  Q.  The  other  difference  is  the  extra  term 
A (7  in  the  RHS  of  Eq.  (2.55).  The  presence  of  this  term  is  due  merely  to  the  fact  that 

one  of  the  equations  is  solved  for  Q  while  the  other  is  solved  for  Q.  This  can  be  seen  by  expanding 
the  term  dQ/dx\ 


3<2  d(J-1)  xdQ 

dx  dx  u  dx  [  ’  dx 

Writing  Eq.  (2.58)  in  discrete  form: 


(2.58) 


Ax  ^  At  1  '  At 


(2.59) 


Multiplying  Eq.  (2.59)  by  At  and  using  a  finite  difference  approximation: 


Q'u  -  Qlk  =  Q"  MJ-1)  +  {J-x)n{Q[k  ~  Q\k )  (2-60) 

Substituting  Eq.  (2.60)  into  Eq.  (2.35),  one  obtains  Eq.  (2.56),  which  is  the  method 
based  on  the  GCL  suggested  by  Thomas  and  Lombard  [49].  Therefore,  Eq.  (2.34)  is  mathematically 
equivalent  to  Eq.  (2.55).  In  order  for  them  to  be  numerically  equivalent,  according  to  the  work  of 
Tamura  and  Fujii  [47],  the  scheme  must  also  be  conservative  in  time.  This  feature  is  achieved  by 
carefully  choosing  the  Jacobian  to  be  used  in  Eq.  (2.35).  They  suggest  using: 
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Figure  2.2.  Schematic  of  the  Spring/Mass/Damper  System 


(^1),  =  5[(^1)"  +  (^1)'1+1]  (2.61) 

B.  STRUCTURAL  DYNAMICS 

Structural  dynamics  modeling  is  based  on  a  two-degree-of-freedom  spring/mass/damper 
system,  shown  in  Fig.  2.2,  to  simulate  the  bending  and  twisting  of  a  wing  section. 

1.  Governing  Equations 

The  equations  governing  a  two-degree-of-freedom  motion  of  structure  are: 


nth  —  Sa  a  +  D/Ji  +  m(ojh  =  L 

(2.62) 

-Sji  +  Ia a  +  Da a  +  Iaora{a  -  a0)  =  M 
where  the  dots  represent  differentiation  with  respect  to  time. 

Equations  (2.62)  are  nondimensionalized  using  reference  length  c,  reference  velocity  doc, 
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reference  mass  poo7t(c/2)2,  and  reference  inertia  pra7t(4~/2)2c2.  Rewriting  Eq.  (2.62)  in  matrix  nota¬ 


tion: 


[M]{X}"  +  [D]{Xy  +  [K}{X}  =  {F} 


(2.63) 


where: 


m  -Sa 

28i,mMookh 

0 

[M}  = 

~Sa  I a 

[D}  = 

0 

2'6(JaM co/^cx 

(2.64) 


{K}  = 

mMlk\  0 

m  =  < 

h 

> 

0  IaMlkl 

a-a0 
.  > 

and 

{F} 

2  , 

=  -Ml{ 

71 

cl 

Cm 

\  / 

> 

where  the  primes  denote  differentiation  with  respect  to  dimensionless  time,  x  =  tam/c,  and  the  other 
parameters  (i.e.,  m,  Ia,  etc.)  are  now  nondimensional. 


2.  Numerical  Technique 

Equation  (2.63)  is  a  system  of  two,  coupled,  second-order,  ordinary  differential  equations. 
Coupling  is  obtained  through  the  terms  containing  Sa  and  the  dependence  of  cl  and  cm  on  h  and 
a.  The  system  is  nonlinear  through  the  nonlinearity  of  cj  and  cm.  Linearization  is  introduced  by 
treating  cl  and  cm  as  constants,  computed  from  the  previous  time-step  of  the  flow  solution. 

Simulations  with  single-degree-of-freedom  may  be  performed  by  setting  Sa  =  0  and  either 
m  =  °°  and  to /,  =  0  or  Ia  =  °°  and  toa  =  0  for  pitching-only  or  plunging-only  motions,  respectively. 
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Equation  (2.63)  is  advanced  in  time  by  solving  for  the  second  derivative  of  the  dependent 


variable  X : 


{X}"  =  [M}-l{F}-[M}-l[K}{X}-[M}-l[D}{Xy  (2.66) 

Rewriting  the  result  as  a  system  of  two,  coupled,  first-order  equations: 

f  {xy  =  {Y} 

<  (2-67) 

( {YY = wrl{n  -  wrl[Km  - 

Finally,  time  integration  is  performed  using  a  first-order  accurate  explicit  Euler  scheme. 
The  use  of  higher  order  methods  does  not  improve  the  quality  of  the  solution  because  the  time 
steps  required  for  the  stability  of  the  unsteady  N-S  equations  yield  very  high  resolution  in  time  for 
Eq.  (2.63). 
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III.  BOUNDARY  CONDITIONS  AND  MULTI-BLOCK 


PROCESSING 

This  chapter  examines  some  background  necessary  for  solving  the  problems  investigated 
in  this  dissertation.  It  is  presented  as  a  separate  chapter  because  some  of  the  material  is  quite 
specific  for  the  layout  of  the  grids  used  throughout  this  investigation.  Furthermore,  this  chapter 
contains  an  important  discussion  about  all  the  boundary  conditions  necessary  for  the  multi-block 
grid  computations.  This  chapter  also  includes  information  that  is  essential  in  order  to  apply  the 
theory  discussed  in  Chapter  II. 

A.  THREE-BLOCK  GRID 

Modeling  of  the  wind  tunnel  geometry  requires  a  mesh  that  is  very  long  in  the  stream  wise 
direction,  compared  to  its  height.  Accurate  computation  of  the  boundary  layer  is  most  easily  per¬ 
formed  on  a  C-grid.  However,  construction  of  a  single -block,  C-grid  would  require  a  highly  skewed 
grid  and  an  excessive  number  of  grid  points  in  regions  far  from  the  airfoil.  Furthermore,  maintaining 
grid  orthogonality  close  to  the  walls  using  a  single -block  grid  would  be  impossible.  These  problems 
can  be  minimized  by  using  a  three-block  grid,  shown  schematically  in  Fig.  3.1,  where  the  governing 
equations  are  solved  in  each  block,  separately. 

A  similar  situation  occurs  with  respect  to  the  airfoil-in-ground-effect  computations.  Be¬ 
cause  the  ground  plane  is  close  to  the  airfoil  section,  the  use  of  a  single -block,  C-grid  would  require 
a  very  skewed  mesh  in  the  vertical  direction.  Once  again,  a  three-block  grid  approach  can  be  used 
instead  of  a  single-block  grid,  as  shown  in  the  schematic  in  Fig.  3.2. 
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Block  (2) 


Block  (3) 


Block  (1) 


Figure  3.1.  Schematic  of  the  Three-Block  Grid  for  the  Flutter  Problem 

The  two  problems  of  this  dissertation  require  two  different  layouts  for  the  three-block  grids. 
The  transonic  flutter  problem  asks  for  the  layout  presented  in  Fig.  3.1.  The  blocks  are  placed  in  line. 
The  first,  or  block  (1),  is  placed  around  the  airfoil.  The  second,  or  block  (2),  discretizes  the  region 
upstream  of  the  first  block  and  the  third,  or  block  (3),  completes  the  downstream  domain.  For  the 
airfoil  in  ground  effect,  the  layout  used  is  shown  in  Fig.  3.2.  The  main  block,  or  block  (1),  represents 
the  region  around  the  airfoil  and  also  the  wake  behind  it.  The  second,  or  block  (2),  fills  the  region 
upstream  of  the  first  one.  The  third,  or  block  (3),  completes  the  domain  above  the  first  two  blocks. 

In  both  problems,  solving  the  N-S  equations  in  all  three  blocks  is  not  necessary.  Blocks  (2) 
and  (3)  are  always  coarse  enough  to  allow  the  use  of  Euler  equations.  Therefore,  N-S  equations  arc 
solved  only  in  block  (1)  for  both  problems. 
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Block  (3) 

Block  (2) 

Block  (1) 

Figure  3.2.  Schematic  of  the  Three-Block  Grid  for  the  Airfoil-in-Ground-Effect  Problem 


B.  BOUNDARY  CONDITIONS 

This  section  describes  the  boundary  conditions  used  for  the  two  different  problems  solved 
in  the  present  work.  The  boundary  condition  for  the  surface  of  the  airfoil  is  the  same  for  both 
problems  and  is  described  next.  The  boundary  conditions  that  are  specific  for  each  problem  arc 
described  in  the  sections  associated  with  each  problem. 

For  inviscid  flow  solutions,  the  viscous  terms  on  the  RHS  of  Eq.  (2.1)  are  set  to  zero,  and  the 
flow-tangency  boundary  condition  is  used  at  the  surface  of  the  airfoil.  For  Navier-Stokes  solutions, 
the  no-slip  condition  is  applied.  Density  and  pressure  arc  extrapolated  to  the  surface  for  both  Euler 
and  Navier-Stokes  solutions. 

For  unsteady  motions,  the  flow-tangency  and  no-slip  conditions  are  modified  to  include 
the  local  motion  of  the  airfoil,  which  also  contributes  to  the  pressure  on  the  surface.  Therefore, 
the  momentum  equation  normal  to  the  surface  (£  direction)  is  solved  to  predict  the  pressure  for  a 
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viscous  flow  more  accurately 


|  I  X  |  wall 

\ p\wall  =  -^  pdf  <  ■  V^  +  ^p|wfl//V^  (3.1) 

y  y  |  waii 

where  x\wau  and  y\waii  are  the  components  of  the  airfoil  velocity.  Furthermore,  Vq  ■  =  0  when 

assuming  that  the  grid  is  orthogonal  at  the  surface.  If  the  airfoil  is  stationary,  the  normal  pressure 
gradient  vanishes  in  agreement  with  boundary-layer  theory. 

1.  Wind-Tunnel  Problem 

The  flow  at  the  wind  tunnel  sections  blocks  (2)  and  (3)  is  assumed  to  be  inviscid  and  the 
Euler  equations  are  solved  in  these  domains.  The  tunnel  walls  are  assumed  to  be  porous,  and 
two  types  of  porous  wall  boundary  conditions  are  applied.  First,  the  porous  boundary  condition  is 
applied  according  to  the  formulation  presented  by  Mokry  et  al.  [35]  (Chapter  2.0  WALL  BOUND¬ 
ARY  CONDITION).  Pressure,  density,  and  x-direction  velocity,  u,  are  extrapolated  from  the  interior 
points  of  the  grid.  The  z-direction  velocity,  w,  is  obtained  according  to  Eq.  (3.2) 

w  P  —  P plenum 

TT  =  a - Tn —  (3-2) 

Uoo  PooU£ 

where  a  is  the  porosity  parameter  of  the  wall  and  p plenum  is  the  pressure  at  the  plenum  chamber.  For 
a  =  0,  for  example,  this  boundary  condition  is  the  same  as  the  flow  tangency  condition  and  the  wall 
is  treated  as  being  completely  solid. 

The  porous-wall  boundary  condition  can  be  further  modified  to  account  for  viscous  effects 
present  at  this  region  of  the  walls  as  follows.  The  tangential  velocity  at  the  wall  is  set  equal  to  zero, 
u  =  0,  and  the  normal  velocity,  w,  is  still  obtained  according  to  Eq.  (3.2).  In  this  work,  this  model 
is  denoted  as  the  porous,  viscous  boundary  condition. 

This  treatment  for  the  porous-wall  boundary  condition  is  quite  different  from  the  one  adopted 
in  the  previous  work  [11].  In  that  study,  the  porous  boundary  was  modeled  by  treating  some  grid 
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Figure  3.3.  Boundary  Conditions  for  the  Flutter  Problem 


cells  as  solid  walls  (applying  flow  tangency)  and  the  others  as  holes  (extrapolating  the  velocity 
vector  from  the  interior  points).  For  the  previous  approach,  the  grid  needed  to  be  approximately 
uniformly  spaced  at  the  porous  region.  Furthermore,  only  certain  discrete  choices  in  the  porosity 
were  allowed.  Values  of  porosity,  such  as  17%,  were  quite  difficult  to  apply  due  to  the  limited 
number  of  grid  cells  along  the  porous  region  of  the  walls.  The  new  model  is  not  dependent  on 
the  streamwise  grid  spacing  at  the  walls  and  also  allows  for  a  continuous  variation  on  the  porosity 
parameter. 

Inflow  and  outflow  boundary  conditions  arc  imposed  in  blocks  (2)  and  (3),  respectively.  For 
the  inflow  boundary,  flow  properties  such  as  pressure,  temperature,  and  velocity  arc  specified  while 
the  density  is  extrapolated  from  the  neighbor  interior  points.  Static  pressure  is  specified  for  the  out¬ 
flow  boundary  condition  and  all  other  properties  are  extrapolated  from  the  interior.  Extrapolations 
in  both  cases  arc  performed  by  using  Riemann  invariants. 

The  boundaries  on  the  left  of  domain  (1)  and  right  of  domain  (2)  are  updated  by  means  of 
an  overlapping  region,  as  shown  in  Fig.  3.3.  There  is  an  overlapping  between  the  blocks  and  the 
flow  variables  arc  transferred  directly  from  one  domain  to  the  other.  The  boundary  conditions  at  the 
right  of  domain  (1)  arc  updated  via  linear-  extrapolation  from  domain  (3)  and  vice-versa. 
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2. 


Airfoil  in  Ground  Effect 


The  flow  at  the  grid  blocks  (2)  and  (3)  is  assumed  to  be  inviscid  and  the  Euler  equations  arc 
solved  in  these  domains.  The  ground  plane  is  assumed  to  be  a  symmetry  plane.  Therefore,  it  can  be 
modeled  as  a  solid  wall  where  the  symmetry  boundary  condition  is  given  by  Eq.  (3.3). 


=  0 

wall 


(3.3) 


where  n  represents  the  normal  direction  at  the  wall. 

The  types  of  boundary  conditions  used  for  all  boundaries  of  the  three-block  grid  are  shown 
in  Fig.  3.4.  Inflow  boundary  conditions  are  applied  to  the  points  on  the  left  boundaries  of  blocks  (2) 
and  (3)  and  also  to  the  upper  boundary  of  block  (3).  Outflow  boundary  conditions  are  specified  to 
the  points  on  the  right  boundary  of  blocks  (1)  and  (3). 

The  boundaries  on  the  left  of  domain  (1)  and  right  of  domain  (2)  are  updated  by  means 
of  an  overlapping  region,  as  shown  in  Fig.  3.4.  There  is  an  overlapping  of  the  blocks  and  the  flow 
variables  are  transferred  directly  from  one  domain  to  the  other.  The  same  type  of  boundary  condition 
is  used  between  blocks  (1)  and  (3),  and  between  (2)  and  (3).  Because  of  the  overlapping  regions, 
direct  transfer  of  the  flow  variables  (p,prr,pw,e)  can  be  used. 


3.  Trailing  Edge  Boundary  Condition 

The  trailing  edge  (TE)  boundary  condition  is  a  crucial  issue  in  both  steady-state  and  un¬ 
steady  computations.  A  schematic  of  a  C-grid  near  the  trailing  edge  is  shown  in  Fig.  3.5.  The  cut 
is  the  region  of  the  computational  domain  corresponding  to  points  for  which  (1  <  i  <  /'/,  k  =  1)  and 
(iu  <  i  <  imax-.k  =1).  These  points  arc  coincident  in  the  physical  domain.  The  trailing  edge  point 
in  the  physical  domain  corresponds  to  the  points  (/  =  ;/)  and  (/  =  iu)  in  the  computational  domain. 
The  points  along  the  cut,  including  the  trailing  edge,  are  shown  separated  in  Fig.  3.5  just  for  clarity. 
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Figure  3.4.  Boundary  Conditions  for  the  Airfoil-in-Ground-Effect  Problem 


TE  -* - cut 


Figure  3.5.  Schematic  of  a  C-Grid  near  the  Trailing  Edge 
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The  boundary  condition  at  the  cut  of  a  C-grid  can  be  done  either  explicitly  or  implicitly. 


In  the  present  method,  this  boundary  condition  is  computed  explicitly,  after  the  flow  tangency  or 
the  no-slip  condition  is  applied  to  the  surface  of  the  airfoil.  The  flow  variables  at  the  cut,  k  —  1, 
arc  calculated  by  averaging  the  upper  and  lower  neighboring  points  for  k  =  2.  Some  options  exist 
when  one  must  decide  how  to  implement  the  explicit  boundary  condition  at  the  trailing  edge.  Three 
possibilities  for  this  implementation  exist: 

•  use  the  averaged  value  from  points  at  (iu. 2)  and  (i/,2); 

•  treat  the  trailing  edge  as  two  different  points  (no  averaging); 

•  use  the  averaged  value  from  points  at  (iu.  1)  and  1). 

The  first  option  is  the  most  obvious  because  it  is  just  an  extension  of  the  boundary  condition 
applied  at  the  cut  of  the  C-grid.  However,  the  first  option  is  not  the  most  suitable  one  because  it 
does  not  enforce  the  no-slip  or  the  flow-tangency  condition  at  the  trailing  edge.  These  boundary 
conditions  arc  overwritten  when  computing  the  average  of  the  off-surface  points  (k  =  2).  In  this 
work,  the  second  alternative  is  called  the  “free  TE  boundary  condition.”  The  no-slip  or  the  flow- 
tangency  condition  is  applied  to  the  TE.  However,  because  the  TE  is  treated  as  two  different  points 
in  the  computational  domain,  a  discontinuity  may  occur  at  the  TE  since  the  two  points  arc  the  same 
in  the  physical  domain.  The  last  option,  called  averaged  TE,  appeal's  to  be  the  most  consistent  one 
for  the  present  work  because,  in  terms  of  discretization  of  the  physical  domain,  the  points  of  the  cut 
corresponding  to  the  TE  are  the  same.  Furthermore,  it  is  guaranteed  that  the  no-slip  condition  or  the 
flow-tangency  condition  is  maintained  at  the  TE  since  these  conditions  are  applied  for  k  =  I  before 
averaging. 

The  influence  of  applying  the  free  or  the  averaged  TE  is  investigated.  A  single  NACA  0012 
airfoil  oscillating  in  a  combined  pitch  and  plunge  motion  is  used  to  determine  the  difference  in  the 
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Figure  3.6.  Grid  around  the  NACA  0012  Airfoil 

solutions  for  both  TE  conditions.  The  reduced  frequency  of  the  motion  is  set  to  k  =  1.0  and  the 
pitch  and  plunge  half  amplitudes  arc  a  =  10  degrees  and  h  =  1.0,  respectively.  Pitch  leads  plunge 
by  (|)  =  90  degrees.  The  pivot  point  for  the  pitch  motion  is  chosen  to  be  the  half  chord  or  xp  =  0.5. 
The  free-stream  Mach  and  Reynolds  numbers  are  =  0.3  and  Re^  =  106,  respectively.  The  B-L 
turbulence  model  is  applied  for  the  fully  turbulent  computations.  The  C-grid  is  presented  in  Fig.  3.6 
and  its  dimensions  are  imax  =  281  and  kmax  =  81. 

The  non-dimensional  equations  of  motion  are  given  by: 

h  =  hsin{Mookx)  a  =  asin(Mockx  +  (}>)  (3.4) 

Steady-state  computations  are  performed  for  a  =  0  degrees,  employing  both  the  averaged 
and  free  TE  boundary  conditions.  The  pressure  distributions  over  the  surface  of  the  airfoil  are 
presented  in  Fig.  3.7.  The  match  is  extremely  good  for  the  two  solutions.  However,  no  conclusion 
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Figure  3.7.  Steady-State  Solutions  for  Averaged  and  Free  TE 

is  drawn  from  the  steady-state  case  because  the  flow  is  symmetric,  and  this  might  be  the  reason  for 
the  excellent  match. 

The  motion  parameters  for  the  oscillation  of  the  airfoil  have  been  chosen  because  they 
force  a  relatively  high  induced  angle  of  attack,  therefore,  allowing  dynamic  stall  to  occur.  The 
time  history  of  the  aerodynamic  coefficients  is  presented  in  Fig.  3.8  for  both  types  of  boundary 
conditions.  Clearly,  the  solution  is  not  periodic.  For  the  same  type  of  boundary  condition,  the  curve 
for  the  lift  coefficients,  for  example,  does  not  behave  periodically. 

The  pressure  distributions  for  various  positions  along  the  cycle  are  shown  in  Fig.  3.9.  The 
left-hand  side  of  this  figure  corresponds  to  positions  along  the  up  stroke  and  the  right-hand  side  to 
the  down  stroke.  Although  some  pressure  distributions  arc  not  in  close  agreement  for  the  two  TE 
conditions,  the  development  of  the  dynamic  stall  is  not  dramatically  different.  In  order  to  better 
clarify  this  point,  entropy  contours  during  the  upstroke  paid  of  the  motion  of  the  airfoil  arc  shown 
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in  Fig.  3.10.  The  vortices  generated  during  the  motion  are  well  captured.  The  left-hand  side  of 


Fig.  3. 10  represents  the  averaged  TE  boundary  condition.  The  free  TE  boundary  condition  is  repre¬ 
sented  on  the  right-hand  side.  Although  the  vortices  do  not  match  perfectly  for  the  two  types  of  TE 
boundary  condition,  their  sizes  are  comparable  and  the  speed  with  which  they  convect  downstream 
is  approximately  the  same.  Note,  also,  that  the  flow  is  not  periodic.  This  can  be  seen  in  Fig.  3.11  in 
which  the  variation  of  the  drag  coefficient  is  shown  against  the  lift  coefficient.  Clearly,  the  solution 
is  not  periodic.  However,  it  stays  inside  a  confined  area  although  it  wanders  from  place  to  place  in 
that  area.  This  behavior  is  represented  in  the  Chaos  Theory  by  the  Torus  Attractor,  which  is  defined, 
for  example,  by  Petree  [38].  This  attractor  appeal's  to  be  the  same  for  both  types  of  TE  boundary 
condition.  Furthermore,  a  computation  of  the  thrust  coefficient  and  propulsive  efficiency  showed 
that  these  parameters  are  virtually  the  same  for  both  TE  conditions  (cj  =  0.40  and  r|  =  0.21).  As 
a  conclusion  from  this  case,  the  type  of  TE  boundary  condition  is  apparently  not  influential  when 
one  is  interested  in  quantities  that  represent  an  average  along  the  cycle.  On  the  other  hand,  the  fine 
details  of  the  flow-field,  as  the  development  of  the  dynamic  stall,  are  strongly  affected  by  the  type 
of  boundary  condition  at  the  trailing  edge. 

A  case  with  no  dynamic  stall  was  also  run  to  evaluate  the  TE  boundary  condition.  A 
NACA  0014  airfoil  oscillating  in  pure  plunge  is  computed.  The  reduced  frequency  of  the  mo¬ 
tion  is  k  =  0.4  with  a  half  amplitude  of  h  =  0.4.  The  free-stream  Mach  and  Reynolds  numbers  are 
M  =  0.3  and  Re  =  106,  respectively.  The  flow  is  assumed  fully  turbulent  and  the  B-L  turbulence 
model  is  used.  The  domain  is  discretized  by  a  321  x  91  C-grid  shown  in  Fig.  3.12. 

When  steady-state  computations  are  performed,  they  deliver  the  same  solution  for  both 
types  of  TE  boundary  condition.  This  is  similar  to  the  NACA  0012  airfoil  case.  Unsteady  simula¬ 
tions  are  also  computed  for  the  two  types  of  TE  boundary  condition.  The  history  of  the  unsteady  lift 
and  drag  aerodynamic  coefficients  is  presented  in  Fig.  3.13.  The  agreement  between  the  solutions 
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Figure  3.10.  Comparison  of  Pressure  Distributions  for  Averaged  and  Free  TE,  NACA  0012,  a  =  10 
degrees,  =  0.3 ,k=  1,  Re <*,  =  1  x  106 
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Figure  3.13.  History  of  Lift  and  Drag  Coefficients  for  the  Pure  Plunge  Case 


for  the  two  different  boundary  conditions  is  quite  good.  The  thrust  coefficient  is  cj  =  0.041  and  the 
propulsive  efficiency  is  r|  =  0.73  for  both  cases.  Therefore,  in  situations  where  there  is  no  dynamic 
stall,  the  difference  of  applying  an  averaged  TE  or  a  free  TE  boundary  condition  is  negligible. 

Although  dynamic  stall  was  computed  somewhat  differently  for  the  two  types  of  TE  bound¬ 
ary  condition,  the  effect  on  the  predicted  thrust  coefficient,  the  propulsive  efficiency,  and  the  con¬ 
vection  speed  of  the  dynamic  vortices  was  minimal.  Therefore,  the  averaged  TE  boundary  condition 
is  used  throughout  the  course  of  this  work. 

Another  important  computation  with  the  averaged  TE  boundary  condition  was  performed  to 
assess  the  influence  of  this  boundary  condition  on  the  dynamic-stall  boundary  suggested  by  Tuncer 
et  al.  [53].  They  computed  the  unsteady  flow  over  a  NACA  0012  airfoil  oscillating  in  pure  plunge 


50 


and  varied  the  reduced  frequency  and  the  amplitude  of  the  motion  to  determine  the  onset  of  dy¬ 
namic  stall.  Their  computations,  using  a  free  TE  boundary  condition,  showed  that  the  dynamic-stall 
boundary  is  given  by  the  equation: 

hk  =  0.35  (3.5) 

Two  computations  with  k  =  1 .0  were  performed  for  h  =  0.3  and  h  =  0.4.  The  first,  cor¬ 
responding  to  hk  =  0.3,  predicted  no  separation  of  the  flow,  in  agreement  with  the  dynamic  stall 
boundary  given  by  Eq.  3.5.  The  second,  corresponding  to  hk  =  0.4,  calculated  a  separated  flow 
during  the  motion  of  the  airfoil,  confirming  that  the  averaged  TE  boundary  condition  also  delivers 
the  dynamic-stall  boundary  described  by  Eq.  3.5. 

A  second  set  of  computations  was  performed  for  the  averaged  TE  boundary  condition  but, 
this  time,  for  the  NACA  0014  airfoil.  The  objective  was  to  determine  the  influence  of  the  shape 
of  the  airfoil  on  the  dynamic-stall  boundary  represented  by  Eq.  3.5.  The  reduced  frequency  was 
kept  equal  to  k  =  1.0  while  the  amplitude  of  the  pure  plunge  motion,  h ,  was  varied.  For  h  =  0.4, 
no  separation  of  the  flow  was  predicted,  suggesting  that  the  dynamic-stall  boundary  might  have 
changed  for  the  NACA  0014  airfoil.  A  computation  with  h  =  0.5  showed  that  the  flow  was  detached 
from  the  airfoil  surface.  These  computations  are  shown  schematically  in  Fig.  3.14,  along  with  the 
dynamic-stall  boundary  of  Eq.  3.5.  Therefore,  clearly,  the  shape  of  the  airfoil  has  a  significant 
influence  on  the  onset  of  dynamic  stall.  For  the  NACA  0014,  the  dynamic-stall  boundary  should  be 
further  investigated  and  it  appeal's  that  it  will  be  given  by: 

hk  >  0.4  (3.6) 
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Reduced  Frequency 


Figure  3.14.  Dynamic-Stall  Boundary  for  the  NACA  0012  Airfoil 

C.  GRID  MOTION 

One  of  the  problems  considered  in  this  dissertation  is  the  numerical  prediction  of  the  flow 
around  an  oscillating  airfoil  inside  a  wind  tunnel.  This  problem  requires  the  use  of  deforming  grids; 
however,  problems  arise  when  adjusting  the  grid,  which  accounts  for  the  motion  of  the  airfoil.  The 
wind  tunnel  walls  are  fixed  at  all  times  while  the  airfoil  is  moving.  Therefore,  the  grid  must  be 
deformed  within  every  time  step  in  order  to  adjust  to  the  relative  movement  between  the  airfoil  and 
the  walls. 

The  required  deformation  of  the  grid  is  obtained  by  dividing  the  whole  domain  of  the  C- 
type  grid  around  the  airfoil  block  (1)  into  four  regions.  These  regions  are  distributed  along  the  main 
block  as  shown  in  Fig.  3.15.  The  first  region  is  called  Ap  and  corresponds  to  the  portion  of  the  block 
that  is  close  to  the  surface  of  the  airfoil  and  that  is  used  to  capture  the  viscous  flow  effects.  In  this 
region,  the  mesh  does  not  deform  but  simply  rotates  and  translates  as  a  solid  body,  following  the 
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same  rotation  and  translation  of  the  airfoil.  This  means  that  no  geometric  change  of  the  grid  cells  in 


the  region  Ap  occurs.  Another  region,  denoted  C,  is  the  one  close  to  the  wind  tunnel  walls  and  to  the 
blocks  (2)  and  (3)  overlapping  sections.  This  region  also  remains  fixed  at  all  times  and,  therefore, 
the  grid  points  do  not  change  in  time  for  an  observer  sitting  on  the  wind  tunnel  walls. 

The  region  Aw  corresponds  to  the  “wake”  following  region  Ap.  The  Aw  is  adjusted  to  the 
movement  of  the  airfoil.  The  adjustment  is  done  using  an  algebraic  grid  generator  that  redistributes 
the  grid  points.  Linear  interpolation  for  the  grid  points  along  a  constant  C,  line  is  applied.  This 
procedure  takes  into  account  that  the  displacement  of  a  grid  point  in  the  region  Aw  is  proportional 
to  the  relative  displacement  of  the  corresponding  points  (same  £  coordinate)  in  regions  A p  and  C. 
Finally,  while  the  relative  location  of  regions  Aw  and  A  p  with  respect  to  the  tunnel  walls  changes, 
region  B  (where  the  grid  deformation  is  the  largest)  is  adjusted  to  provide  a  smooth  grid  between 
regions  A  and  C.  This  adjustment,  again,  is  accomplished  by  the  algebraic  grid  generator,  but  now 
lineal-  interpolation  is  performed  along  constant  £,  lines. 

The  same  approach  is  applied  for  block  (1)  of  the  airfoil-in-ground-effect  configuration. 
Blocks  (2)  and  (3)  stay  fixed  throughout  the  unsteady  computations  and  block  (1)  must  be  adjusted 
to  the  relative  movement  between  the  airfoil  and  the  ground  plane.  This  adjustment  is  done  exactly 
as  described  for  the  wind-tunnel  case. 
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Figure  3.15.  Schematic  of  the  Regions  for  Grid  Motion 

D.  PARALLEL  PROCESSING 

The  strategy  of  dividing  the  domain  into  blocks  also  has  the  advantage  of  allowing  the  use 
of  parallel  processing.  Each  one  of  the  three  blocks  is  solved  in  a  separate  computer  and  data  from 
the  grid  points  located  on  the  common  boundaries  arc  exchanged. 

Data  related  to  the  boundaries  of  each  block  must  be  transferred  to  the  neighboring  blocks 
to  be  used  as  boundary  conditions.  Data  are  transferred  by  a  standard  library  called  Message  Passing 
Interface  (MPI).  The  standard  library  used  in  the  present  work  is  the  Local  Access  Multimachine 
(MPI/LAM)  from  the  University  of  Notre  Dame. 

The  parallel  version  of  the  Navier-Stokes  solver  used  in  the  present  work  is  based  on  the 
code  validated  and  used  by  several  other  researchers  (for  instance:  Ekaterinaris  and  Menter  [16], 
Ekaterinaris  and  Platzer  [17],  Sanz  and  Platzer  [43],  and  Tuncer  and  Platzer  [51]).  Then,  subroutines 
dedicated  only  for  data  passing  between  neighboring  blocks  were  added  to  the  previous  version  of 
the  code  to  make  it  parallel  processing  capable. 
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The  blocks  are  defined  during  the  grid-generation  phase.  They  do  not  need  to  be  of  equal 


dimensions.  Actually,  because  N-S  equations  are  solved  in  block  (1)  and  Euler  equations  are  solved 
in  blocks  (2)  and  (3),  the  first  block  tends  to  be  larger  than  the  other  two.  This  characteristic  of 
the  parallel  code  allows  the  use  of  different  machines  to  compose  the  cluster.  For  example,  block 
(2)  is  usually  the  smallest  one.  If  block  (1)  is  an  N-S  block  with  300x100  points  and  block  (2)  is 
60x50,  the  former  is  10  times  larger  than  the  latter.  Therefore,  one  can  use  a  machine  much  slower 
to  compute  block  (2)  as  the  machine  used  for  block  (1). 

A  task  performed  in  the  present  work  was  evaluating  the  ability  of  the  parallel  version  of 
the  solver  to  reduce  the  computational  time  required  for  convergence  of  a  particular  problem.  A 
simulation  for  an  NLR  7301  airfoil  inside  a  solid- wall  wind  tunnel  was  run  and  a  reduction  of 
approximately  10%  in  the  time  required  to  achieve  5,000  iterations  was  found.  Note  that  the  three- 
block  grid  used  in  this  case  was  not  favorable  for  a  substantial  reduction  in  the  computational  time, 
since  the  size  of  the  first  block  corresponded  to  approximately  85%  of  the  total  number  of  grid 
cells.  Furthermore,  thin  layer  Navier-Stokes  equations  were  solved  in  it  while  Euler  equations  were 
solved  in  the  other  two  blocks.  Consequently,  a  10%  reduction  in  the  computational  time  is  quite 
satisfactory  for  this  particular  case  and  demonstrates  the  feasibility  of  the  parallel  version  of  the 
code.  The  pressure  distribution  obtained  with  the  parallel  code  for  a  steady-state  problem  is  shown 
in  Fig.  3.16.  This  case  is  a  NFR  7301  airfoil  inside  a  wind  tunnel  with  solid  walls.  It  can  be  seen 
that  the  parallel  solution  matches  the  single  processor  solution  well. 

An  unsteady  case  was  also  run  in  order  to  check  the  parallel  version  of  the  code.  The  results 
of  unsteady  computations  for  a  porous  wind  tunnel  wall  case  are  shown  in  Fig.  3. 17.  The  agreement 
between  the  parallel  and  single  processor  solutions  is  also  very  good. 

For  the  parallel  version  of  the  code,  the  three  blocks  are  solved  simultaneously,  each  one 
in  a  different  node  of  the  PC  cluster.  The  values  of  the  flow  variables  at  the  boundaries  arc  always 
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Figure  3.17.  Comparison  of  History  of  Angle  of  Attack  for  Porous  Walls 
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exchanged  among  the  nodes  before  the  beginning  of  each  iteration,  for  all  blocks.  Therefore,  the 


exchanged  boundary  conditions  are  all  at  the  same  computational  time.  On  the  other  hand,  for  the 
single  processor  solution,  the  blocks  are  solved  sequentially.  The  boundary  conditions  for  block  (2), 
for  example,  arc  transferred  from  block  (1)  after  its  solution  is  advanced  one  time  step.  This  means 
that  the  boundary  conditions  are  not  exactly  at  the  same  computational  time.  Despite  this  subtle 
difference  in  the  way  the  two  codes  work,  they  yielded  virtually  the  same  solution.  This  shows  that, 
for  a  reasonably  small  time  step,  whether  the  boundary  conditions  are  exchanged  at  the  exact  same 
computational  time  or  not  is  not  influential. 
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IV.  VALIDATION  OF  THE  MODELS 


A.  DEFORMING  GRIDS 

The  problems  studied  in  this  work  need  a  solver  that  can  handle  deforming  grids.  The  first 
question  that  naturally  comes  to  mind  is  how  this  methodology  handles  constant  changes  in  the 
coordinates  of  the  grid  points.  In  order  to  answer  this  question,  two  tests  arc  designed.  They  entail 
a  subsonic  and  a  transonic  flow,  respectively.  Both  tests  show  that  the  present  method  can  perform 
computations  with  deforming  grids  with  a  negligible  deviation  of  the  solution. 

The  next  two  subsections  explain  in  detail  how  these  two  tests  are  performed  and  the  results 
obtained  in  each  test  arc  shown. 

1.  Subsonic  Flow  Test 

This  test  is  designed  to  demonstrate  the  ability  of  this  method  to  perform  unsteady  com¬ 
putations  with  a  deforming  grid  for  subsonic  flow.  The  problem  is  the  pure  plunge  oscillation  of  a 
NACA  0014  airfoil.  The  airfoil  is  moving  with  a  free-stream  Mach  number  =  0.3  and  a  reduced 
frequency  k  =  0.5.  The  half  amplitude  of  the  plunge  motion  is  0.4c. 

First,  an  unsteady  inviscid  solution  is  performed  using  a  single  grid  around  the  airfoil.  This 
grid  moves  along  with  the  airfoil  as  a  solid  body.  Therefore,  the  grid  is  moving  but  not  deforming. 
The  grid  used  for  the  computations  is  a  C-grid  with  281  points  in  the  c,  direction  and  51  in  the  C, 
direction.  The  single-block  grid  is  shown  in  Fig.  4.1. 

Second,  another  solution  is  computed  using  a  four-block  grid.  The  main  block  is  a  C-grid 
that  contains  the  airfoil.  The  second  block  is  an  H-grid  placed  just  upstream  of  the  main  block.  The 
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Figure  4.1.  Close-Up  of  the  Single-Block  Grid  around  the  Airfoil 


third  block  is  placed  right  on  top  of  the  first  two.  The  fourth  block  is  a  mirror  image  of  the  third 
with  respect  to  the  chord  plane  of  the  airfoil.  This  grid  is  shown  in  Fig.  4.2.  Blocks  two,  three,  and 
four  stay  fixed  at  all  times.  The  only  block  allowed  to  deform,  in  order  to  accommodate  the  motion 
of  the  airfoil,  is  the  main  one.  The  outer  boundary  of  this  block  stays  fixed  in  order  to  exchange 
boundary  conditions  with  the  other  blocks. 

The  comparison  of  the  two  solutions  is  shown  in  Fig.  4.3.  The  solution  for  the  single-block 
grid  is  represented  by  the  dashed  line  and  the  solution  for  the  four-block  grid  corresponds  to  the  solid 
line.  Despite  a  very  large  deformation  of  the  first  block  in  the  four-block  solution,  the  agreement 
with  the  single-block  solution  is  quite  reasonable.  The  maximum  deviations  for  the  lift,  drag,  and 
moment  coefficients  arc  around  2.1%,  3.5%,  and  2.7%,  respectively. 

The  upper  and  lower  boundaries  of  the  first  block  for  the  four-block  grid  were  placed  at 
h  =  0.7  and  li  =  —0.7,  respectively.  The  highest  and  lowest  positions  for  the  airfoil  correspond  to 
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Figure  4.2.  Portion  of  the  Four-Block  Grid  around  the  Airfoil 


h  =  0.4  and  h  =  —0.4,  respectively.  This  means  that  the  deformation  of  the  first  block  is  quite  large. 
To  better  illustrate  this  point,  the  four-block  grid  with  the  airfoil  in  its  lowest  position  ( h  =  —0.4) 
is  shown  in  Fig.  4.4.  A  steady-state  solution  for  the  airfoil  at  h  =  —0.4  is  evaluated  and  compared 
to  a  similar  computation  for  h  =  0  (undeformed  grid)  in  Fig.  4.5.  The  solver  is  not  able  to  predict 
a  symmetric  solution  for  h  =  —0.4,  as  it  is  for  h  =  0,  due  to  the  large  deformation  of  the  grid. 
This  deviation  of  the  steady-state  solution  for  h  =  —0.4  may  explain  the  difference  between  the 
single-block  and  the  four-block  solutions  for  the  unsteady  problem. 

Despite  the  large  deformation  of  the  grid,  this  method  was  able  to  keep  the  difference  be¬ 
tween  the  single -block  and  the  four-block  solutions  lower  than  4%.  It  is  important  to  note  that  the 
four-block  grid  used  here  is  just  an  extension  of  the  three-block  grid  used  to  solve  the  airfoil-in- 
ground-effect  problem. 
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Figure  4.6.  Three-Block  Grid  near  the  NLR  7301  Airfoil 


2.  Transonic  Flow  Test 

This  test  is  intended  to  determine  the  performance  of  this  method  for  a  transonic  flow  prob¬ 
lem  in  which  only  the  grid  is  changing.  Two  different  solutions  are  computed.  One  for  a  fixed  grid 
and  another  for  a  deforming  grid. 

The  problem  solved  is  the  transonic  flow  around  the  NLR  7301  airfoil  inside  a  wind  tunnel. 
The  flow  is  steady  at  a  Mach  number  of  M  =  0.768,  a  Reynolds  number  of  Re  =  1.727  x  106,  and  an 
angle  of  attack  of  a  =  1 .28  degrees.  The  walls  of  the  wind  tunnel  arc  considered  as  being  completely 
solid.  The  grid  used  is  shown  in  Fig.  4.6.  It  is  a  three-block  grid  and  the  dimensions  of  the  first, 
second,  and  third  blocks  are  281  x  81,  41  x  41,  and  41x61,  respectively. 

To  validate  the  handling  of  the  deforming  grid,  the  test  is  performed  where  the  boundary 
of  region  B  in  block  (1)  remains  fixed,  but  the  interior  points  of  B  perform  sinusoidal  oscillations. 
Grid  points  in  blocks  (2)  and  (3)  as  well  as  the  ones  in  regions  Ap,  Aw,  and  C  of  block  (1)  remain 
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Figure  4.7.  Deviation  from  the  Fixed-Grid  Solution 


fixed.  The  reduced  frequency  of  the  oscillation  of  grid  points  in  region  B  is  set  to  k  =  0.3  and  the 
maximum  amplitude  to  Ah  =  0.05. 

This  test  is  initialized  from  a  steady-state  solution  for  a  solid  wind  tunnel  wall  boundary 
condition  at  a  non-dimensional  time  of  x  =  46.72.  Ideally,  the  solution  should  remain  constant 
after  the  oscillation  of  the  grid  points  in  region  B  has  started.  Nonetheless,  because  of  numerical 
errors  and  the  movement  of  grid  points  across  a  shock  wave,  some  variation  of  the  aerodynamic 
coefficients  should  be  expected.  The  results  for  the  test  arc  shown  in  Fig.  4.7.  The  amplitudes  of 
the  deviation  of  the  aerodynamic  coefficients  arc  Ac[  =  3.2  x  10  4,  Ac  d  =  1.6  x  1 0  5 ,  and  A cm  = 
6.4  x  10-5. 

The  results  of  a  DFT  analysis  on  the  deviations  is  presented  in  Fig.  4.8.  One  can  see  that 
the  non-dimensional  frequency  of  the  signal  is  computed  as  /  =  0.0375.  This  corresponds  to  a 
reduced  frequency  of  k  =  0.307,  which  is  approximately  the  reduced  frequency  in  which  the  grid 
was  oscillating.  The  amplitudes  of  the  oscillations  are  approximately  three  orders  of  magnitude 
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Figure  4.8.  DFT  Analysis  of  the  Deviation 

smaller  than  the  oscillations  of  the  aerodynamic  coefficients  during  the  limit-cycle  oscillations. 
These  small-amplitude  oscillations  are  considered  to  be  acceptable. 

This  result  for  the  error  between  the  fixed  and  deforming  grid  solution  was  found  to  be 
consistent  with  the  results  obtained  in  similar  tests  performed  in  [47] .  This  means  that  the  present 
method  is  capable  of  handling  deforming  grids  with  an  acceptable  degree  of  accuracy. 

B.  PRESCRIBED  MOTION  ANALYSIS 

First  in  this  section,  the  results  of  the  present  solver  arc  compared  with  computations  of 
Isogai  et  al.  [21]  and  Tuncer  et  al.  [53]  for  a  single  NACA  0012  in  a  pitch/plunge  motion.  The  pa¬ 
rameters  of  the  motion  of  the  airfoil  were  chosen  in  such  a  way  that  separation  occurs  and  dynamic 
vortices  are  developed.  The  second  paid  of  this  section  contains  a  comparison  with  the  unsteady  po¬ 
tential  flow  solver  USPOT  for  an  NACA  0014  airfoil  in  ground  effect.  Euler  solutions  are  computed 
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and  the  good  agreement  with  predictions  of  USPOT  shows  that  the  present  code  can  handle  moving 


and  deforming  grids  with  a  reasonable  degree  of  accuracy.  The  results  also  show  a  significant  influ¬ 
ence  of  the  free-stream  Mach  number  on  the  aerodynamic  coefficcients  computed  for  the  airfoil  in 
ground  effect. 

1.  Single-Airfoil  Pitch/Plunge  Motion 

The  results  of  the  present  code  arc  compared  to  the  ones  obtained  by  Isogai  et  al.  [21]  for  a 
pitch/plunge  motion  of  a  NACA  0012  airfoil.  The  reduced  frequency  is  constant  and  equal  to  k  =  1.0 
for  all  computations.  The  half  amplitudes  arc  a  =  10  degrees  and  h—  1.0  for  the  pitch  and  plunge 
motions,  respectively.  The  phase  angle  between  pitch  and  plunge  motions  is  varied  from  tf>  =  30 
degrees  to  (f)  =  150  degrees.  The  free-stream  Mach  and  Reynolds  numbers  are  set  to  —  0.3  and 
Rea o  =  105,  respectively.  The  fully  turbulent  flow  is  computed  by  using  the  B-L  turbulence  model. 
The  grid  used  is  the  same  presented  in  Fig.  3.6. 

The  vorticity  contours  are  presented  in  Fig.  4.9  for  (|)  =  90  degrees.  This  case  differs  only 
by  the  free-stream  Reynolds  number  from  the  one  used  in  Section  III.B.3,  in  which  the  trailing  edge 
boundary  condition  is  discussed.  The  results  obtained  by  Isogai  et  al.  [21]  are  on  the  LHS  and 
the  present  results  arc  on  the  RHS  of  Fig.  4.9.  The  agreement  between  the  two  solutions  is  good. 
The  location  of  the  primary  vortex  is  predicted  very  closely  in  both  solutions.  The  prediction  of  a 
secondary  vortex  occurs  in  both  cases.  Even  a  small  separation  near  the  trailing  edge  is  picked  up 
by  both  solvers.  Note  that  the  contours  are  not  the  same  for  both  results.  The  contour  lines  of  Isogai 
et  al.  are  not  defined  in  [21],  hence,  it  is  difficult  to  reproduce  the  same  contours. 

In  their  research,  Isogai  et  al.  [2 1]  compared  their  results  with  Tuncer  et  al.  [53] .  They  varied 
the  phase  angle  between  pitch  and  plunge  and  calculated  the  thrust  coefficient  and  the  propulsive 
efficiency.  Their  results  for  these  two  parameters  are  compared  to  the  present  code  computations  in 
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Figure  4.9.  Vorticity  Contours  of  Isogai  et  al.  (left)  and  Present  Work  (right) 


Fig.  4.10.  In  the  present  work,  the  reduced  frequency  of  the  motion  is  equal  to  k  =  1.0,  the  plunge 
amplitude  is  h  =  1.0,  and  the  Reynolds  numbers  are  Re „  =  105  and  Re <*>  =  106.  Note  that,  because 
of  a  different  nondimensionalization,  the  reduced  frequency  and  the  plunge  amplitude  of  [21]  arc 
k  =  0.5  and  h  =  2,  respectively. 

The  present  results,  shown  in  Fig.  4.10,  contain  error  bars  due  to  the  non-periodicity  of  the 
computed  solution.  The  thrust  coefficient  and  the  propulsive  efficiency  are  not  the  same  for  each 
cycle  and  the  variation  of  these  parameters  through  ten  cycles  is  represented  by  the  error  bars.  It  is 
important  to  note  that  the  maximum  induced  angle  of  attack  is  symmetric  with  respect  to  a  phase 
angle  equal  to  (|)  =  90  degrees.  For  instance,  the  maximum  induced  angle  of  attack  for  (|)  =  30 
degrees  is  the  same  as  for  (|)  =  150  degrees.  The  present  results  show  a  similar  behavior  for  the 
thrust  coefficient  and  propulsive  efficiency.  The  curve  through  the  mean  values  of  these  para  maters 
is  almost  symmetric  with  respect  to  cf)  =  90  degrees,  although  the  variation  around  the  mean  value 
is  not  symmetric. 

An  attempt  to  compare  results  with  Ramamurti  and  Sandberg  [40]  was  made  but,  because 
of  the  large  values  of  amplitudes  of  the  motion  and  reduced  frequency,  the  present  code  was  unable 
to  compute  a  meaningful  solution.  The  value  of  the  Strouhal  number  used  in  their  computation  is 
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Figure  4.10.  Thrust  Coefficient  and  Propulsive  Efficiency  for  the  NACA  0012  Airfoil 
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around  Sr  =  0.6,  which  corresponds  to  a  reduced  frequency  k  =  3.8.  The  amplitudes  of  the  motion 
are  a  =  15  degrees  and  h  =  1.0.  With  this  combination  of  parameters,  the  present  code  predicts 
locally  supersonic  flow  even  for  a  free-stream  Mach  number  of  AW  =  0.1.  Because  Ramamurti  and 
Sandberg  [40]  used  an  incompressible  Navier-Stokes  solver  for  their  computations,  no  comparison 
with  their  results  could  be  done. 

2.  Airfoil-in-Ground-Effect  and  Pure-Plunge  Motion 

A  final  validation  of  the  present  code  is  performed  for  a  biplane  configuration.  This  compu¬ 
tation  shows  that  the  equations  for  deforming  grids  can  compute  the  flow  over  an  oscillating  airfoil 
in  ground  effect  accurately. 

Euler  solutions  for  a  pure  plunging  NACA  0014  biplane  arc  compared  with  the  one  ob¬ 
tained  by  using  an  unsteady  potential  flow  solver  called  USPOT  [37],  The  reduced  frequency  of 
the  oscillation  is  k  =  0.5.  The  half  amplitude  of  the  motion  is  h  =  0.4.  The  separation  between  the 
wings  is  /  =  1.4.  For  the  present  solution,  only  one  wing  is  modeled  and  the  presence  of  the  second 
wing  is  simulated  by  a  symmetry  plane.  The  distance  from  the  wing  to  the  symmetry  plane  is,  then, 
d  =  0.7. 


The  Euler  solutions  were  run  for  Mach  numbers  of  AW  =  0.2  and  AW  =  0.3.  The  multi¬ 
block  grid  used  is  shown  in  Fig.  4.11.  The  dimensions  for  blocks  (1),  (2),  and  (3)  are  289  x  41, 
41  x  39,  and  165  x  51,  respectively.  The  solutions  for  the  two  Mach  numbers  are  presented  in 
Fig.  4.12  and  compared  with  the  USPOT  solution.  The  incompressible  flow  solution  computed  by 
USPOT  corresponds  to  a  Mach  number  of  AW  =  0.  The  Euler  solutions  for  AW  =  0.2  and  AW  =  0.3 
show  a  maximum  lift  coefficient  higher  than  the  one  predicted  by  USPOT.  Nonetheless,  clearly,  the 
tendency  is  to  approach  the  incompressible  solution  as  the  Mach  number  diminishes.  The  same 
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Figure  4.11.  Euler  Multi-Block  Grid  near  the  Airfoil 


behavior  is  observed  with  respect  to  the  drag  coefficient.  The  minimum  values  for  the  Euler  drag 
coefficients  tend  to  move  toward  the  USPOT  solution  as  the  Mach  number  decreases. 

The  influence  of  the  Mach  number  on  the  aerodynamic  coefficients  is  known  to  be  quite 
significant  for  wings  in  extreme  ground  effect,  as  discussed  by  Rozhdestvensky  [41].  For  wings 
in  extreme  ground  effect,  meaning  distances  from  the  ground  lower  than  10%  of  the  chord,  the 
aerodynamic  coeffients  vary  with  the  parameter  (1  — Moo)-1  and  not  (1  —Moo)-1/2.  Even  though  the 
smallest  distance  from  the  airfoil  to  the  ground  during  the  oscillation  is  cl  =  0.3.  the  influence  of  the 
Mach  number  is  still  significant. 

This  computation  shows  that  the  solution  obtained  with  the  compressible,  deforming-grid 
code  is  in  good  agreement  with  the  solution  calculated  for  the  unsteady  potential  flow.  Therefore, 
the  present  solver  shows  that  it  can  perform  accurately  for  a  deforming  grid  situation.  Furthermore, 
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Figure  4. 12.  Comparison  of  Euler  and  Potential  Flow  Solutions  for  a  Biplane  Configuration 


the  solution  shows  a  significant  influence  of  the  Mach  number  on  the  unsteady  flow  characteristics 
of  the  biplane  or  airfoil-in-ground-effect  configuration. 
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V.  RESULTS 


A.  FLAPPING  AIRFOIL  PROPULSION 

The  biplane  model  studied  at  the  Naval  Postgraduate  School  is  presented  in  Fig.  5.1.  It 
has  two  sets  of  wings  with  a  chord  of  c  =  0.064 m.  The  wings  oscillate  sinusoidally  in  an  opposed 
plunge  manner.  The  mean  separation  between  the  two  planes  is  /  =  0.0896m,  which  yields  a  non- 
dimensional  value  of  l  —  l/c  —  1.4.  This  configuration,  as  explained  in  Chapter  III,  can  be  modeled 
by  only  one  airfoil  oscillating  above  a  symmetry  plane,  similar  to  a  wing  in  ground  effect.  The 
non-dimensional  mean  distance  from  the  airfoil  section  to  the  symmetry  plane  is,  then,  d  =  0.7. 
The  amplitude  of  the  motion  of  the  airfoil  corresponds  to  a  non-dimensional  value  of  h  =  0.4.  This 
enables  a  minimum  clearance  of  d  —  h  =  0.3  from  the  symmetry  plane. 

This  configuration  was  tested  by  Lund  [34]  and  the  experiments  were  conducted  inside 
a  low-speed,  open-section  wind  tunnel  with  speeds  varying  from  £/«,  =  0  to  £/«,  =  9.5m/ s.  The 
Reynolds  number  varied  from  Re „  =  0  to  Rerx,  =  45 . 000  and  the  Mach  number  from  M„  =  0  to  M.x  = 
0.028.  The  frequency  of  the  oscillation  of  the  wings  was  set  to  discrete  values  of  /  =  3, 5,  and  7  Hz. 

The  numerical  computations  were  performed  in  a  three-block  grid  in  which  block  (1)  corre¬ 
sponds  to  a  289  x  81  C-grid,  block  (2)  is  a  61  x  29  H-grid,  and  block  (3)  is  represented  by  a  152  x  51 
H-grid.  A  close-up  of  the  grid  around  the  airfoil  is  presented  in  Fig.  5.2.  The  S-A  and  B-L  turbu¬ 
lence  models  were  used  to  perform  fully  turbulent  calculations.  Fully  laminar  computations  were 
also  performed  because  of  the  low  Reynolds  number  of  the  experiments.  In  this  case,  no  turbulence 
model  was  used. 
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notch  in  rear  nacelle 


Figure  5.2.  Three-Block  Grid  near  the  Airfoil  Section 
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No  sensitivity  analysis  to  the  size  of  the  grid  was  performed  in  the  present  work  because  of 


the  results  obtained  by  Tuncer  et  al.  [53].  They  found  that  the  computed  solution  was  not  sensitive  to 
grid  sizes  of  121  x  62,  241x61,241x91,  and  311  x  71.  However,  a  sensitivity  study  to  the  time  step 
was  performed  for  some  of  the  unsteady  cases  of  the  present  work.  This  study  was  accomplished 
by  increasing  the  number  of  Newton  sub-iterations  within  each  physical  time  step,  which  was  found 
to  be  equivalent  to  reducing  the  time  step.  The  results  of  this  procedure  showed  that  the  parameters 
computed  by  means  of  averaging  through  a  cycle  did  not  change  significantly.  The  prediction  of  fine 
detail  of  the  flow-field,  such  as  the  frequency  of  leading  edge  vortex  shedding  and  the  location  of  the 
dynamic  vortex  along  the  cycle,  was  more  sensitive  to  the  time  step.  Nevertheless,  because  the  main 
interest  of  the  present  work  is  on  averaged  values  along  the  cycle,  such  as  the  thrust  coefficient  and 
the  propulsive  efficiency,  three  Newton  sub-iterations  were  used  and  the  time  step  was  kept  within 
a  minimum  of  3,500  steps  per  cycle. 

1.  Steady-State  Computations 

Steady-state  solutions  arc  needed  to  start  up  unsteady  solutions.  Another  reason  for  com¬ 
puting  them  is  the  need  to  correct  the  calculation  of  the  thrust  coefficient.  In  the  experiments,  one 
is  interested  in  measuring  the  thrust  due  to  the  oscillation  of  the  wings.  This  is  done  by  measur¬ 
ing  the  drag,  or  thrust,  for  the  entire  model,  including  the  wings  and  fuselage,  and  subtracting  the 
steady-state  drag  measured  for  the  same  model.  The  difference  is  the  thrust  due  to  the  flapping  of 
the  wings. 

The  consequence  of  this  procedure  is  that  all  types  of  steady-state  drag  arc  eliminated, 
including  friction  and  pressure  drag  of  the  wing.  Hence,  in  order  to  compare  numerical  and  experi¬ 
mental  results,  the  same  procedure  must  be  performed.  The  computed  steady-state  drag  also  has  to 
be  subtracted  from  the  unsteady  values. 


75 


Table  5.1.  Steady-State  Drag  Coefficients 


h 

laminar 

S-A 

S-A 

B-L 

M  =  0.1 

M  =  0.1 

§ 

II 

p 

u> 

M  —  0.1 

-0.4 

0.0260 

0.0038 

0.0038 

0.0050 

-0.2 

0.0228 

0.0030 

0.0030 

0.0031 

0.0 

0.0217 

0.0028 

0.0028 

0.0029 

0.2 

0.0212 

0.0027 

0.0027 

0.0028 

0.4 

0.0209 

0.0026 

0.0026 

0.0027 

For  the  biplane  configuration,  or  wing  in  ground  effect,  it  is  known  that  the  drag  coefficient 
is  a  function  of  the  ground  clearance.  Consequently,  a  series  of  steady-state  computations  were 
performed  for  various  values  of  distances  from  the  ground.  The  results  for  these  calculations  is 
presented  in  Table  5.1. 

The  steady-state  solutions  for  fully  turbulent  flows  were  run  for  a  Reynolds  number  Rerx,  = 
106.  This  value  is  much  higher  than  the  Reynolds  number  of  the  experiment  but,  in  order  to  use  B-L 
or  S-A  turbulence  models,  106  is  the  minimum  value  of  their  range  of  application.  The  pressure 
distributions  for  these  turbulence  models  arc  shown  in  Fig.  5.3  for  a  distance  from  the  symmetry 
plane  equals  to  d  =  0.3,  corresponding  to  the  closest  position  during  the  sinusoidal  motion  of  the 
airfoil.  The  two  solutions  arc  very  close  to  each  other,  showing  that,  for  steady-state  computations, 
the  turbulence  models  do  not  have  a  strong  influence  on  the  calculations.  Furthermore,  the  pressure 
distribution  is  not  symmetric  anymore.  The  lift  force  resulting  from  the  imbalance  of  pressure  is 
now  in  the  direction  of  the  ground  plane. 
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Figure  5.3.  Pressure  Coefficient  Distributions  for  B-L  and  S-A  Models 


The  solution  for  fully  laminar  flow  was  performed  using  Re „  =  104.  The  pressure  distribu¬ 
tion  over  the  airfoil  for  this  case  is  presented  in  Fig.  5.4.  Also  shown  in  this  figure  is  the  pressure 
distribution  for  the  fully  turbulent  flow  using  the  S-A  turbulence  model.  The  two  solutions  look 
quite  different.  The  flow  for  the  laminar  solution  is  separated  whereas  the  fully  turbulent  flow  is 
always  attached  to  the  airfoil.  A  snapshot  of  the  vorticity  field  for  the  laminar  computation  is  pre¬ 
sented  in  Fig.  5.5.  The  boundary  layer  is  quite  thick  and  detaches  from  the  surface  of  the  airfoil, 
producing  vortices  that  arc  shed  from  the  trailing  edge.  This  means  that  even  the  fixed-airfoil  so¬ 
lution  has  an  unsteady  behavior  for  the  low  Reynolds  number  in  question.  This  unsteady  behavior 
can  also  be  seen  in  the  history  of  aerodynamic  coefficients  shown  in  Fig.  5.6. 

The  vortices  produced  by  the  detachment  of  the  boundary  layer  arc  drag  inducing  vortices. 
Their  orientation  is  such  that  they  tend  to  approach  the  ground  plane.  Therefore,  the  wake  is  de¬ 
flected  downward  and  induces  a  positive  (upward)  lift  on  the  airfoil. 
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Figure  5.6.  History  of  Aerodynamic  Coefficients  for  the  Laminar  Flow  Solution 


The  steady-state  solutions  show  that  the  fully  turbulent  flow  is  not  influenced  significantly 
by  the  use  of  the  B-L  or  the  S-A  turbulence  model.  The  flow  is  always  attached  to  the  airfoil  and 
the  lift  force  is  negative  (pointing  to  the  ground  plane).  On  the  other  hand,  the  laminar-flow  solution 
yields  a  relatively  thick  boundary  layer,  which  detaches  from  the  airfoil  and  produces  an  oscillatory 
behavior  for  the  aerodynamic  coefficients. 

2.  Oscillating  Airfoil  Computations 

In  the  experiments,  the  frequency  of  oscillation  of  the  wings  was  set  to  a  fixed  value,  and 
the  speed  of  the  wind  tunnel  was  changed  in  order  to  vary  the  reduced  frequency.  This  means 
that  the  free-stream  Reynolds  number  is  different  for  each  value  of  reduced  frequency.  The  range 
of  the  free-stream  Reynolds  number  in  the  experiments  was  from  Re„  =  0  to  Re „  —  45,000.  In 
the  numerical  solutions,  the  Reynolds  number  was  kept  fixed  for  all  values  of  reduced  frequency 
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because  of  the  large  number  of  steady-state  solutions  that  would  have  to  be  generated  to  start  up  the 
unsteady  computations.  The  Reynolds  number  was  set  equal  to  Re <*,  =  104  for  the  laminar  solutions 
and  Rea 0  =  106  for  the  fully  turbulent  computations.  It  is  important  to  mention  that  this  procedure 
may  induce  discrepancies  when  comparing  experimental  and  numerical  data. 

The  free-stream  Mach  number  is  another  important  issue.  Ideally,  because  of  the  limitations 
of  the  numerical  scheme,  the  minimum  free-stream  Mach  number  should  be  around  =  0.3. 
However,  for  this  Mach  number  and  values  of  reduced  frequency  around  k  =  2.0,  the  local  Mach 
numbers  of  the  flow-field  reach  values  as  high  as  M  =  1.5,  at  which  point  compressibility  effects 
arc  severe.  The  Mach  contour  lines  for  these  conditions  arc  shown  in  Fig.  5.7  with  the  airfoil  in 
a  position  very  close  to  h  =  0  during  the  up  stroke.  The  contour  lines  go  from  M  =  0  to  M  =  1.6 
with  an  increment  of  AM  =  0.04.  The  region  where  the  local  flow  is  supersonic  is  presented  in 
Fig.  5.8(a).  Although  this  region  is  not  very  large,  the  region  where  compressibility  effects  arc 
already  significant  is  much  larger.  The  entropy  contour  lines  are  shown  in  Fig.  5.8(b).  It  can  be  seen 
that  the  supersonic  pocket  occurs  inside  the  dynamic  vortex  which  is  being  generated  at  the  leading 
edge.  Hence,  the  physics  of  the  dynamic  stall  are  completely  altered  by  compressibility  effects.  It 
is  also  important  to  mention  that  the  TVD  scheme  is  switched  off  for  the  airfoil  in  ground-effect 
computations. 

A  computation  for  a  free-stream  Mach  of  =  0.1  was  performed  for  the  same  reduced 
frequency  of  k  =  2.0.  The  Mach  contour  lines  are  presented  in  Fig.  5.9  with  the  same  increment 
of  AM  =  0.04.  The  local  Mach  number  for  this  case  stays  below  M  =  0.55  for  the  whole  cycle. 
Therefore,  most  of  the  computations  were  performed  for  a  free-stream  Mach  number  of  =  0.1. 
The  idea  was  to  eliminate  any  source  of  compressibility  effects  since  the  Mach  numbers  of  the 
experiments  were  much  lower. 
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Figure  5.9.  Mach  Contour  Lines  for  =  0.1  and  k  =  2.0 

The  computed  thrust  coefficient  for  the  airfoil  in  ground  effect  configuration  and  Ma »  =  0.1 
is  presented  in  Fig.  5.10  as  a  function  of  the  reduced  frequency.  The  turbulence  models  used  for  the 
fully  turbulent  solutions  were  B-L  and  S-A.  The  Reynolds  number  for  these  cases  was  Re„  =  106. 
For  the  fully  laminar  solution,  the  Reynolds  number  was  Re^  =  104.  The  potential  flow  solutions  for 
a  biplane  (USPOT)  and  a  single  airfoil  (UPOT)  configurations  and  the  experimental  values  obtained 
by  Lund  [34]  in  the  NPS  wind  tunnel  are  also  shown  in  the  same  figure. 

The  thrust  coefficient  per  airfoil  is  defined  by  Eq.  (5.1): 

CT  =  TJl^Ulc  (5‘1} 

where  T  =  —D  is  the  thrust  force  generated  by  the  flapping  motion  of  the  airfoil,  is  the  free- 
stream  density  and  LL,  is  the  free-stream  velocity  of  the  flow. 

Computations  for  =  0.3  were  also  performed  for  the  S-A  turbulence  model.  They  arc 
compared  with  the  solutions  for  the  same  model  and  =  0. 1  in  Fig.  5.11.  For  =  0.3  and  values 
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Figure  5.10.  Thrust  Coefficient  for  =  0.1 


of  reduced  frequency  higher  than  k  =  1.0,  compressibility  effects  clearly  change  the  behavior  of  the 
solution  dramatically. 

The  agreement  with  the  experiments  of  [34]  is  much  better  for  the  laminar  than  for  the  fully 
turbulent  solutions,  regardless  of  the  turbulence  model  used  for  turbulent  flow.  This  is  expected 
because  the  Reynolds  numbers  of  the  experiments  are  much  closer  to  the  Reynolds  number  of  the 
laminar  solution  than  to  the  Reynolds  number  of  the  fully  turbulent  computations.  At  this  range 
of  Reynolds  numbers,  the  vortices,  which  arc  shed  due  to  the  detachment  of  the  boundary  layer, 
dominate  the  flow  over  the  airfoil,  especially  for  low  frequencies.  In  fact,  solutions  for  S-A  and  B-L 
turbulence  models,  when  no  vortices  arc  shed,  agree  better  with  the  USPOT  solution.  This  means 
that  the  prediction  of  the  vortices  produced  by  boundary  layer  detachment  or  the  prediction  of  the 
dynamic  stall  is  crucial  for  simulating  correctly  the  behavior  of  the  flow  at  this  range  of  Reynolds 
numbers. 
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Figure  5.11.  Thrust  Coefficient  for  the  S-A  Turbulence  Model 


The  entropy  distribution  of  the  laminar  flow  field  around  the  airfoil  is  presented  in  Fig.  5.12 
through  a  series  of  frames  for  different  positions  along  the  cycle.  The  LHS  of  Fig.  5.12  represents 
the  down  stroke  whereas  the  RHS  shows  the  up  stroke.  The  dynamic  stall  vortices  generated  during 
the  motion  arc  well  captured  by  the  solver,  and  it  is  clear  that  the  boundary  layer  is  quite  thick  for 
some  positions  along  the  cycle.  Therefore,  reviewing  the  thin  layer  assumption  of  the  present  solver 
is  essential.  This  assumption  may  not  hold  for  positions  along  the  cycle  if  the  boundary  layer  is 
relatively  thick. 

The  solutions  for  turbulent  flow  arc  close  to  the  potential  flow  computations  (USPOT)  for 
reduced  frequencies  up  to  k  =  1.0  (S-A  and  M0 <,  =  0.3)  and  k  =  2.0  (S-A,  B-L,  and  M0 0  =  0.1). 
As  stated  by  Tuncer  et  al.  [53],  the  limit  for  separated  flows  for  a  single  NACA  0012  airfoil  is 
given  by  the  equation  hk  =  0.35.  According  to  the  discussion  presented  in  Section  III.B.3,  this  limit 
could  be  further  extended  to  hk  =  0.45  for  the  NACA  0014  airfoil.  For  h  =  0.4,  the  value  of  the 
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reduced  frequency  for  flow  separation  becomes  k  =  1.125.  However,  the  computed  thrust  coefficient 


remains  close  to  the  potential  flow  solution  even  beyond  this  dynamic  stall  limit.  In  order  to  assess 
this  behavior,  two  sequences  of  entropy  contours  for  k  =  2.0  and  M„  =  0.1  arc  given  in  Figs.  5.13 
and  5.14  for  the  S-A  and  B-L  turbulence  models,  respectively.  For  the  S-A  solution.  Fig  5.13, 
the  flow  is  not  completely  separated.  A  recirculation  bubble,  which  convects  downstream  but  only 
separates  from  the  surface  of  the  airfoil  near  the  trailing  edge,  is  predicted.  Furthermore,  the  bubble 
predicted  for  the  upper  surface  is  larger  than  the  bubble  computed  for  the  lower  surface.  For  the 
B-L  simulation.  Fig.  5.14,  the  dynamic  stall  vortex  is  clearly  predicted  by  the  solver  for  both  upper 
and  lower  surfaces.  These  solutions  suggest  that  the  turbulence  model  influences  the  prediction  of 
dynamic  stall,  and  also  suggest  that  the  dynamic  stall  limit  for  the  biplane  configuration  may  be 
slightly  different  from  the  single-wing  value. 

The  propulsive  efficiency  is  defined  as  the  power  generated  by  the  thrust  force  divided  by  the 
power  necessary  to  oscillate  the  airfoil,  and  it  is  represented  in  Eq.  (5.2).  The  computed  propulsive 
efficiency  is  shown  in  Fig.  5.15. 


fl  = 


(5.2) 


where  T  =  —D  is  the  thrust  force  and  Preq  is  the  power  required  for  oscillating  the  airfoil. 

It  is  known  that  the  biplane  configuration  produces  more  thrust  per  wing  than  the  single 
wing.  A  comparison  between  computed  values  of  propulsive  efficiency  and  thrust  per  wing  for  the 
biplane  and  single  configurations  is  shown  in  Fig.  5.16.  It  can  be  seen  that  the  airfoil  in  ground 
effect  produces  more  thrust  per  wing,  with  almost  the  same  efficiency,  than  the  single  airfoil.  This 
behavior  was  demonstrated  previously  by  Jones  et  al.  [24] . 
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Figure  5.14.  Entropy  Contour  Lines  for  B-L  Turb.  Model,  =  0. 1,  and  k  =  2.0 
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Figure  5.15.  Propulsive  Efficiency  for  the  Biplane  Configuration 
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(a)  (b) 

Figure  5.16.  Comparison  between  the  Biplane  and  the  Single -Airfoil  Configurations 


B.  TRANSONIC  FLUTTER 


A  major  conclusion  by  Castro  et  al.  [11]  was  identifying  how  significant  the  porosity  of 
the  wind-tunnel  walls  was  on  predicting  transonic  steady  and  unsteady  flow  characteristics  of  the 
airfoil  inside  the  wind  tunnel.  The  current  work  was  conducted  in  order  to  further  investigate  the 
effect  of  porosity  of  the  walls  and  to  improve  the  method  of  prescribing  such  a  boundary  condition. 
Therefore,  the  approach  adopted  was  based  on  determining  a  porous  boundary  condition  that  yielded 
a  steady-state  surface  pressure  distribution  for  the  NLR  7301  airfoil,  in  the  presence  of  the  wind- 
tunnel  walls,  which  best  agreed  with  the  experiments  of  Schewe  et  al.  [31].  Once  a  satisfactory 
porous  boundary  condition  was  determined,  the  flutter  computations  were  performed.  Adjusting 
other  flow  parameters  of  the  experiment,  such  as  inflow  and  outflow  boundary  conditions,  was  not 
attempted  since  they  were  not  given  in  [31].  The  wind-tunnel  test  was  performed  at  a  Mach  number 
of  Moo  =  0.768,  an  angle  of  attack  of  a  =  1.28  degrees,  and  a  Reynolds  number  of  Re  =  1.727  x  106. 
The  same  flow  conditions  were  used  in  the  present  simulation,  since  no  corrections  were  applied  to 
the  experimental  data  [31].  The  plenum  pressure  was  always  kept  equal  to  the  free-stream  pressure 
(p plenum  =  P<=»).  The  results  for  the  steady  calculations  were  used  as  the  start  for  the  unsteady 
simulations. 

All  steady-state  and  unsteady  computations  were  performed  using  a  C-type  281  x  81  point 
main-block  grid  (1),  shown  in  Fig.  5.17,  which  was  generated  from  the  original  NLR  7301  air¬ 
foil  surface  data.  Blocks  (2)  and  (3)  were  Cartesian-type  and  contained  41  x  41  and  41  x  61 
points  along  the  streamwise  and  the  normal  directions,  respectively.  The  Spalart-Allmaras  (S-A) 
and  Baldwin-Lomax  (B-L)  turbulence  models  were  used  throughout  the  course  of  the  unsteady 
computations  for  the  NLR  7301  airfoil. 

All  the  computations  were  performed  in  a  time-accurate  mode  using  a  constant  time  step. 
For  fixed  angles  of  incidence,  the  solution  was  run  for  a  long  time  after  convergence  in  residuals 
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Figure  5.17.  C-Type  Grid  near  the  NLR  7301  Airfoil 


so  that  variations  in  all  aerodynamic  coefficients  were  reduced  to  machine  zero.  This  was  done  in 
order  to  ensure  that  all  flow  disturbances  from  the  initial  transients  were  swept  out  of  the  domain. 
At  convergence,  all  solutions  at  fixed  angles  of  incidence  did  not  exhibit  any  unsteadiness. 

A  grid-sensitivity  study  was  performed  previously  by  Weber  et  al.  [55]  showing  that  results 
were  not  significantly  changed  for  grid  sizes  of  221  x  91  with  an  initial  wall  spacing  of  2  x  10  5, 
which  yields  y+  <  1.0,  with  40  points  in  the  wake  and  with  the  far- field  boundary  placed  20  chord 
lengths  away  from  the  airfoil  surface.  Therefore,  an  initial  wall  spacing  of  1  x  10~5  was  chosen  to 
keep  y+  ss  1.0  even  for  unsteady  computations.  The  grid  size  of  281  x  81  chosen  for  block  (1)  of 
the  present  computations  guarantees  that  there  will  be  40  points  in  the  wake.  Furthermore,  because 
the  far- field  boundary  is  placed  at  the  tunnel  walls  and  not  20  chord  lengths  away  from  the  airfoil, 
8 1  grid  points  in  the  ^  direction  are  sufficient  to  guarantee  a  reasonable  grid  resolution  away  from 
the  airfoil. 
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A  time-sensitivity  analysis  was  conducted  in  the  present  work.  The  time  step  used  in  the 


present  investigation  was  such  that  a  minimum  of  2,400  iterations  per  cycle  of  oscillation  was 
applied  to  each  one  of  the  unsteady  cases.  The  time  step  of  some  unsteady  computations  was 
divided  by  two,  and  no  significant  change  in  the  flutter  frequency  or  inter-modal  phase  angle  was 
observed.  However,  a  reduction  of  approximately  5%  on  the  flutter  amplitudes  occurred  for  these 
cases.  As  seen  later  in  this  section,  other  parameters  influenced  the  flutter  amplitudes  much  more 
strongly  than  the  time  step.  Therefore,  most  of  the  effort  was  concentrated  in  studying  the  influence 
of  these  parameters  rather  than  the  influence  of  the  time  step.  Furthermore,  the  computational  cost 
of  using  a  finer  time  resolution  for  most  of  the  unsteady  simulations  would  be  enormous. 

1.  Steady-State  Computations 

First,  a  computation  for  a  solid  wind  tunnel  wall  was  performed  using  the  S-A  turbulence 
model.  The  normal  component  of  the  velocity  at  the  walls  was  set  to  zero  and  the  other  flow 
variables  were  extrapolated  from  the  interior  points  of  the  grid.  The  pressure  distribution  for  this 
case  is  presented  in  Fig.  3.16  and  compared  with  the  averaged  pressure  distribution  of  Schewe  et 
al.  [31].  The  agreement  with  the  experiment  can  be  further  improved,  as  discussed  in  the  previous 
work  [1 1],  by  modeling  the  porosity  of  the  wind-tunnel  walls. 

The  model  for  the  porosity  of  the  wall,  based  on  the  theory  presented  by  Mokry  et  al.  [35], 
Eq.  (3.2),  was  evaluated  next.  The  pressure  distributions  for  25%  and  50%  wall  porosity  parameters 
are  presented  in  Fig.  5.18.  The  porous,  inviscid  boundary  condition  was  used  for  these  cases.  It  can 
be  seen  that  a  better  agreement  with  the  experiment  is  achieved  for  a  porous  boundary  condition 
than  for  the  solid  wall.  That  the  turbulence  model  has  a  significant  influence  on  the  prediction  of 
the  shock  location  is  also  clear. 
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Figure  5.18.  Surface  Pressure  Distribution  for  a  Porous,  Inviscid  Boundary  Condition  at  the  Tunnel 
Walls 
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Figure  5.19.  Surface  Pressure  Distribution  for  a  Porous,  Viscous  Boundary  Condition  at  the  Tunnel 
Walls 


The  porous,  viscous  model  was  also  tested  for  porosity  parameters  of  a  =  0.25  and  a  =  0.50 
and  yielded  slightly  better  agreement  with  the  experiment  than  the  porous,  inviscid  model.  The 
comparison  of  the  pressure  distributions  with  the  experimental  results  for  this  case  are  shown  in 
Fig.  5.19.  The  porosity  parameter  of  a  =  0.25  seems  to  yield  a  better  agreement  with  the  time- 
averaged  pressure  distribution  of  [31],  predicting  a  better  location  for  the  upper  shock.  Once  again, 
the  turbulence  model  significantly  affects  the  computed  shock  location. 

In  summary,  the  steady-state  computations  demonstrate  that  the  porosity  of  the  walls  has 
a  very  strong  influence  on  the  computed  flow  held  around  the  NLR  7301  airfoil.  Additionally,  the 
method  of  applying  the  corresponding  boundary  condition  and  the  turbulence  model  are  influential. 
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2. 


Flutter  Computations 


Based  on  the  conclusions  drawn  from  the  steady  computations,  the  flutter  simulations  were 
performed  using  the  porous  model  for  the  wind-tunnel  walls.  The  nominal  conditions  of  the  wind- 
tunnel  test  were  preserved,  namely  =  0.768,  a,-  =  1.28  degrees,  and  Re  =  1.727  x  106,  as  well  as 
the  spring-neutral  angle  of  attack  ao  =1.91  degrees.  Note  that  these  values  were  modified  in  [55], 
to  account  for  wind-tunnel  interference,  to  =  0.753,  a ,-c  =  —0.08  degrees.  Re  =  1.727  x  106, 
and  aoc  =  0.635  degrees.  All  time-accurate  flutter  computations  were  performed  assuming  fully 
turbulent  flow  using  the  Spalart-Allmaras  and  Baldwin-Lomax  turbulence  models.  The  computa¬ 
tions  predicted  flutter  in  two  degrees  of  freedom.  Limit-cycle  oscillations  (LCO)  were  computed  in 
agreement  with  the  wind-tunnel  test. 

In  the  experimental  test  case  [31],  limit-cycle  oscillations  in  pitch  and  plunge  were  reported. 
The  experiment  was  conducted  at  a  total  pressure  of  0.45  bar  and  a  dynamic  pressure  of  0.126  bar. 
A  time-averaged  angle  of  attack  of  a  =  1.28  degrees  was  measured  for  an  angle  of  attack  at  wind- 
off  condition  of  ao  =  1.91  degrees,  which  is  equivalent  to  the  spring-neutral  angle  of  attack  in 
the  numerical  simulation.  The  porosity  parameter  associated  with  the  perforated  wall  of  the  DLR- 
Gottingen  wind  tunnel  was  stated  as  a  =  0.25.  The  holes  in  the  wall  were  drilled  at  an  angle  of  30 
degrees  with  respect  to  the  surface  of  the  wall.  No  measurements  of  the  pressure  at  the  wall  were 
performed;  therefore,  the  plenum  pressure  was  assumed  to  be  the  free-stream  pressure  in  the  present 
study.  The  dimensionless  structural  parameters  of  the  experiment  are  summarized  in  Table  5.2.  The 
same  parameters  were  used  for  the  aeroelastic  computation. 

Initially,  the  flutter  computations  were  started  based  on  a  porosity  parameter  of  a  =  0.25. 
Time  histories  of  the  angle  of  attack  for  a  porous,  inviscid  boundary  condition  at  the  tunnel  walls 
are  shown  in  Figs.  5.20  and  5.21  for  the  S-A  and  B-L  turbulence  models,  respectively.  Note  that 
LCO  has  not  been  achieved  for  both  turbulence  models.  The  initial  oscillations  damp  out  and  the 
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Table  5.2.  Structural  Parameters 


xp  =  0.2500 

ka  =  0.3330 

xa  =  0.0484 

kh  =  0.2540 

m  =  946.00 

8a  =  0.0041 

4  =  33.900 

8*  =  0.0073 
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Figure  5.20.  Angle-of-Attack  History  for  a  =  0.25,  Porous,  Inviscid  BC,  and  S-A  Turb.  Model 


computations  converge  to  a  steady  value  of  the  angle  of  attack.  For  the  B-L  model,  the  oscillations 
arc  erratic  after  some  time,  showing  a  different  behavior  from  the  S-A  turbulence  model. 

Next,  a  porous,  viscous  boundary  condition  was  evaluated.  Time  histories  of  the  angle 
of  attack  for  this  case  arc  presented  in  Fig.  5.22  and  5.23.  For  a  porous,  viscous  model,  LCO  is 
clearly  achieved  for  both  the  S-A  and  B-L  turbulence  models,  but  the  amplitudes  are  higher  than 
the  reported  experimental  values  of  Knipfer  et  al  [31].  The  flutter  results  for  the  computations  arc 
summarized  in  Table  5.3.  Differently  from  the  porous,  inviscid  model,  the  behavior  of  the  solutions 
for  the  porous,  viscous  boundary  condition  is  similar  for  both  turbulence  models.  No  erratic  oscil- 
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Figure  5.21.  Angle-of-Attack  History  for  a  =  0.25,  Porous,  Inviscid  BC,  and  B-L  Turb.  Model 


lations  are  present  and  the  flutter  parameters  are  predicted  closely.  Therefore,  no  significant  effect 
of  the  turbulence  model  on  computations  with  the  porous,  viscous  model  existed. 

The  computations  for  a  =  0.25  show  that  the  type  of  boundary  condition  for  the  porous 
wall  can  lead  to  different  unsteady  solutions.  LCO  is  obtained  for  the  porous,  viscous  boundary 
condition  in  which  the  flow  is  essentially  vertical  at  the  wall,  but  for  the  porous,  inviscid  boundary 
condition,  the  flow  is  almost  tangent  to  the  wall  and  the  motion  is  damped.  In  fact,  neither  the  vis¬ 
cous  nor  the  inviscid  boundary  condition  accurately  models  the  details  of  the  near-wall  flow  of  the 
experiment.  The  experiment  used  holes  drilled  at  30  degrees  imposing  a  curvature  to  the  flow.  The 
approach  adopted  in  the  present  work  uses  two  extreme  conditions  of  flow  curvature  represented 
by  the  porous,  inviscid  and  viscous  boundary  conditions  (tangent  and  normal,  respectively).  The 
curvature  of  the  flow  in  the  experiment  is  between  these  two  extreme  values.  Fortunately,  the  exper¬ 
imental  LCO  amplitudes  also  lie  between  the  computed  amplitudes  for  the  two  types  of  boundary 
condition.  This  fact  suggests  that  the  porous  boundary  condition  for  the  wind-tunnel  walls  sig¬ 
nificantly  improves  the  computations  of  LCO  for  the  NLR  7301  inside  the  DLR-Gottingen  wind 
tunnel. 
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Table  5.3.  Flutter  Results 


method 

a 

a 

h 

/ 

[deg] 

[deg] 

[mm] 

[Hz] 

[deg] 

Exp  .a 

1.28 

0.18 

0.65 

32.85 

176.7 

S-A  b 

1.19 

0.00 

0.00 

34.4 

166 

S-Ac 

1.15 

1.70 

4.68 

34.5 

165 

B-L  i, 

1.11 

0.00 

0.00 

- 

- 

B-Lc 

0.98 

1.79 

5.00 

34.6 

165 

S-Arf 

1.24 

0.78 

2.9 

36.7 

149 

S-A, 

0.07 

3.78 

11.1 

32.30 

171.8 

a  =  without  wind-tunnel  corrections. 
b  =  with  porous,  inviscid  wall;  a  =  0.25. 
c  =  with  porous,  viscous  wall;  a  =  0.25. 


d  =  previous  work;  with  50%  porosity. 
e  =  [55];  fully  turbulent  (unbounded  computation). 
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Figure  5.24.  Mach  Contour  Lines  for  Porosity  a  =  0.25  and  Porous,  Inviscid  Boundary  Condition 

The  Mach  contour  lines  of  the  flow  field  near  the  airfoil  for  the  porous,  inviscid  and  porous, 
viscous  boundary  conditions  arc  presented  in  Figs.  5.24  and  5.25,  respectively.  The  Mach  contours 
lines  go  from  M  =  0  to  M  =  2.0  with  an  increment  of  AM  =  0.05  for  both  cases.  The  flow  goes 
to  a  relatively  low  speed  near  the  porous  walls  for  the  porous,  viscous  boundary  condition  case, 
as  shown  in  Fig.  5.25.  For  the  porous,  inviscid  boundary  condition,  the  flow  near  the  tunnel  walls 
remains  at  a  relatively  high  speed,  as  shown  in  Fig.  5.24. 

In  order  to  study  the  influence  of  the  porosity  parameter  of  the  wind-tunnel  walls  in  pre¬ 
dicting  limit-cycle  oscillations,  a  parametric  variation  of  a  was  conducted.  First,  a  porous,  inviscid 
boundary  condition  was  used  to  compute  a  solution  for  a  porosity  parameter  of  50%.  The  time  his¬ 
tory  of  the  angle  of  attack  shows  that  the  LCO  was  never  achieved  and  the  initial  oscillations  damp 
out  as  illustrated  in  Fig.  5.26.  Next,  several  cases  were  run  for  different  values  of  the  porosity  pa¬ 
rameter  using  the  porous,  viscous  boundary  condition.  The  LCO  amplitudes  obtained  in  those  runs 
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Figure  5.26.  Angle-of-Attack  History  for  a  =  0.50,  Porous,  Inviscid  BC,  and  S-A  Turb.  Model 


are  shown  in  Fig.  5.27.  These  results  are  also  provided  in  Table  5.4.  That  the  porosity  parameter 
also  has  a  significant  influence  on  the  LCO  amplitudes  is  clear.  Because  of  the  flow  curvature  at  the 
porous  walls,  the  LCO  amplitudes  are  expected  to  be  smaller  than  the  ones  presented  in  Fig.  5.27 
and  Table  5.4.  Nevertheless,  the  influence  of  the  porosity  parameter  is  believed  to  be  the  same. 
Final  results  for  the  flutter  computations,  including  flutter  frequency,  phase,  amplitudes  a  and  h , 
and  mean  angle  of  attack  a,  are  shown  in  Tables  5.3  and  5.4. 

The  influence  of  the  tunnel  blockage  was  also  investigated.  Three  more  cases  were  run  with 
different  heights  of  the  wind-tunnel  test  section,  H.  The  nominal  condition  for  the  experiment  was 
H  =  1.0 m,  which  yields  H jc  =  3.33.  The  additional  cases  represented  H  jc  —  5.00,  H  jc  =  6.67, 
and  H/c  =  °°  (unbounded  flow).  All  nominal  values  of  the  experiment  were  preserved.  Note  that  the 
unbounded  solution,  in  this  case,  is  different  from  the  results  obtained  in  [55],  due  to  the  corrected 
free-stream  conditions  they  used.  The  history  of  the  angle-of-attack  amplitude  for  these  cases  can 
be  seen  in  Fig.  5.28.  The  porous,  viscous  boundary  condition  with  a  =  0.25  was  used  at  the  tunnel 
walls  in  all  the  cases  except  the  unbounded  flow.  The  results  show  a  tendency  of  decreasing  the 
LCO  amplitudes  as  the  height  of  the  test  section  H  is  increased.  For  H  jc  =  5.00,  the  oscillations 
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Figure  5.28.  Variation  of  LCO  Amplitudes  with  the  Solid  Blockage 


are  already  damped  out  as  in  the  unbounded  flow  solution.  As  expected,  the  results  tend  to  approach 
the  unbounded  flow  solution  as  H /c  is  increased.  These  results  also  show  that,  depending  on  the 
value  of  H /c  chosen  for  the  wind-tunnel  test,  the  flutter  characteristics  of  the  measurements  may 
differ  significantly  from  the  free-flight  situation.  Moreover,  apparently,  the  variation  of  the  LCO 
amplitudes  with  the  tunnel  blockage  is  rather  non-linear. 

The  influence  of  the  Mach  number  on  the  flutter  characteristics  of  the  NLR  7301  airfoil 
including  tunnel  walls  was  studied  as  well.  In  order  to  save  computational  efforts,  the  porous, 
inviscid  wall  boundary  condition  was  used.  This  type  of  wall  boundary  condition  requires  a  much 
shorter  computational  time  to  achieve  LCO.  The  time  history  of  the  angle-of-attack  amplitudes 
for  some  Mach  numbers  is  presented  in  Fig.  5.29  for  both  tunnel  and  unbounded  solutions.  The 
amplitudes  for  wind-tunnel  solutions  were  always  higher  than  the  amplitudes  for  unbounded  flow 
in  the  Mach  number  range  of  the  computations. 

The  results  obtained  in  this  work  show  that  the  porosity  parameter,  the  solid  blockage  of 
the  test  section,  and  the  Mach  number  all  have  a  strong  influence  on  the  prediction  of  the  transonic 
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Figure  5.29.  Variation  of  Angle-of- Attack  Amplitudes  with  Mach  Number 


flutter  characteristics  of  the  NLR  7301  airfoil,  as  indicated  in  Figs.  5.20  through  5.29.  The  type 


of  boundary  condition  used  for  the  porous  wind-tunnel  wall  is  also  influential.  The  computed  am¬ 
plitudes  of  a  and  h  are  larger  than  the  ones  measured  in  the  experiment  [31]  for  a  porous,  viscous 
boundary  condition  but  lower  for  a  porous,  inviscid  boundary  condition.  The  type  of  boundary  con¬ 
dition  has  only  a  small  effect  on  predicting  the  flutter  frequency  and  the  inter-modal  phase  angle. 
These  values  were  predicted  more  closely  in  the  present  work  than  in  [11].  Therefore,  it  is  con¬ 
sidered  that  an  improvement  was  obtained  with  respect  to  the  previous  work  [11].  The  predicted 
flutter  frequency  deviates  from  the  experimental  value  by  4.9%  and  the  inter-modal  phase  angle 
by  1 1  degrees.  Nonetheless,  these  parameters  were  predicted  more  closely  by  the  unbounded  flow 
computations  of  Weber  et  al  [55].  The  values  of  frequencies,  amplitudes,  and  phase  angles  were 
calculated  by  means  of  a  DFT-analysis  of  the  last  10  cycles. 

Although  the  amplitudes  were  overpredicted  or  underpredicted  depending  on  the  type  of 
wall  boundary  condition,  the  limit-cycle  oscillation  phenomenon  was  correctly  predicted  and  the 
frequency  and  the  inter-modal  phase  angle  were  computed  within  a  reasonable  accuracy.  It  must 
be  kept  in  mind  that  the  plenum  chamber  pressure  was  assumed  to  be  the  free-stream  pressure. 
Also,  inevitably,  other  uncertain  factors  exist.  It  appeal's,  then,  that  the  wind-tunnel  porosity  model 
used  in  the  present  computations  significantly  improves  the  prediction  of  flutter  characteristics  of 
wind-tunnel  flows. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A  parallel  version  of  the  solver  was  implemented  and  compared  with  the  single-processor 
version.  The  parallel  code  reduced  the  wall-clock  time  required  for  a  solution  and  yielded  identical 
results  to  the  single-processor  solution  for  both  steady-state  and  unsteady  computations.  A  cluster 
of  Linux  PC’s  was  used  to  conduct  the  computations  for  both  problems  under  investigation  in  this 
work. 

The  reason  for  studying  the  airfoil-in-ground-effect  problem  is  to  understand  the  physics  of 
the  low-speed  flow  over  the  opposed-plunge  biplane  of  the  NPS’s  micro-air  vehicle.  When  using 
a  compressible  solver  for  high  reduced  frequencies,  the  local  Mach  number  may  reach  supersonic 
values  at  which  compressibility  effects  are  dominant,  as  discussed  in  Chapter  V.  Since  the  micro-air 
vehicle  will  fly  at  much  lower  Mach  regimes,  predictions  using  compressible  solvers  may  not  cap¬ 
ture  the  correct  physics  of  this  problem.  Furthermore,  as  indicated  by  Anderson  [3],  it  appeal's  that 
dynamic-stall-vortex  capture  at  high  reduced-frequencies  makes  thrust  generation  more  efficient. 
Hence,  in  order  to  investigate  numerically  motions  with  high  reduced-frequencies,  the  use  of  an 
incompressible  solver  would  be  more  appropriate  because  it  would  avoid  compressibility  effects. 

Another  important  issue  is  the  thin-boundary-layer  assumption.  Micro-air  vehicles  fly  at 
low  Reynolds  numbers  and,  in  such  case,  this  assumption  is  questionable  because  boundary  layers 
are  relatively  thick  at  this  regime.  Therefore,  it  is  more  appropriate  to  use  an  incompressible,  full 
Navier-Stokes  code  to  study  numerically  the  unsteady  flow  over  the  biplane  wing  of  the  NPS’s 
micro-air  vehicle  for  low  Mach  and  Reynolds  numbers. 

The  Spalart-Allmaras  turbulence  model  was  modified  so  that  it  would  work  with  moving 
and  deforming  grids.  This  change  is  important  because  this  model  was  developed  for  fixed  grids 
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and  the  convection  of  vorticity  was  not  correctly  predicted  for  moving  or  changing  grids.  This 


change  alone  is  insufficient  to  work  with  moving  and  deforming  grids.  The  turbulence  model  must 
also  be  solved  with  the  same  time  step  as  the  conservation  equations.  In  some  solvers,  a  relaxation 
scheme  is  used  in  which  the  time  step  for  the  turbulence  model  is  different  from  the  time  step  used 
for  the  conservation  equations.  Since  relaxation  is  used  in  the  present  solver,  the  Spalart-Allmaras 
turbulence  model  must  be  further  modified  to  be  time  accurate. 

It  was  found  that  assuming  fully  laminar  flow  for  the  airfoil-in-ground-effect  problem  de¬ 
livered  the  best  agreement  with  the  experiments  conducted  by  Lund  [34].  Despite  the  limitations  of 
the  solver,  the  agreement  of  the  numerical  predictions  with  the  experiments  was  good  for  reduced 
frequencies  up  to  k  =  2.0.  Beyond  this  value,  compressibility  effects  played  a  major  role  in  the  nu¬ 
merical  predictions.  Therefore,  computations  were  not  attempted  for  values  of  reduced  frequency 
beyond  k  =  2.0. 

A  study  of  the  trailing-edge  boundary  condition  revealed  that  the  manner  in  which  it  is 
implemented  in  the  flow  solver  is  influential  when  there  are  non-linearities  in  the  flow  field.  For 
unsteady  flows  with  low  induced  angles  of  attack,  where  the  flow  is  always  attached  to  the  airfoil, 
almost  no  difference  was  found  in  the  computations  when  using  the  free  or  the  averaged  TE  bound¬ 
ary  conditions.  However,  if  the  induced  angle  of  attack  was  high  enough  to  generate  dynamic  stall, 
the  computed  solution  was  sensitive  to  the  TE  boundary  condition.  When  non-linearities  are  present, 
the  predicted  solution  is  no  longer  periodic.  Instead,  the  predicted  solution  follows  an  attractor,  as 
shown  in  Chapter  III.  Although  the  fine  detail  of  the  solution  is  dependent  on  the  type  of  TE  bound¬ 
ary  condition,  the  attractor  followed  by  the  solution  is  apparently  not  significantly  influenced.  A 
similar  behavior  was  found  regarding  the  influence  of  the  turbulence  model.  When  non-linearities 
arc  present,  the  solution  is  significantly  affected  by  the  use  of  the  S-A  or  the  B-L  turbulence  models. 
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In  the  present  solver,  the  boundary  conditions  are  computed  explicitly.  The  grid-cut  bound¬ 


ary  conditions  arc  obtained  by  means  of  linear  interpolation  and  the  airfoil  surface  boundary  con¬ 
ditions  are  performed  using  a  first-order  extrapolation  for  Euler  computations  and  zeroth-order  ex¬ 
trapolations  for  Navier-Stokes  calculations.  An  improvement  on  the  accuracy  of  the  solver  could  be 
obtained  by  computing  the  cut  boundary  conditions  implicitly  and  using  higher-order  extrapolations 
for  the  airfoil  surface  boundary  conditions. 

An  improved  method  of  modeling  the  porosity  of  the  wind  tunnel  walls  was  implemented. 
Two  different  models  were  tested.  One  adopted  the  approach  presented  by  Mokry  et  al.  [35],  in 
which  the  normal  velocity  through  the  porous  region  was  proportional  to  the  pressure  difference  be¬ 
tween  the  plenum  chamber  and  the  test  section.  The  other  consisted  of  a  viscous  approach  assuming 
the  same  normal  velocity  but  no  tangential  component  at  the  porous  region  of  the  wind-tunnel  walls. 
The  modeling  of  the  tunnel  wall  porosity  and  the  way  the  corresponding  boundary  condition  was 
applied  (viscous  or  inviscid)  were  both  found  to  affect  the  numerical  predictions  of  the  steady-state 
and  flutter  characteristics  significantly.  The  improved  porosity  models  also  allowed  more  flexibil¬ 
ity  for  the  generation  of  the  grid  because  the  requirement  for  an  almost  equally  spaced  grid  at  the 
porous  region  of  the  walls  was  no  longer  necessary. 

The  transonic  two-degree-of-freedom  bending/torsion  flutter  analysis  of  the  NLR  7301  su¬ 
percritical  airfoil  section  was  performed  with  tunnel  walls  modeled  with  an  improved  porosity 
boundary  condition.  This  model  showed  that  the  porosity  parameter  influences  significantly  the 
prediction  of  the  limit-cycle  amplitudes.  On  the  other  hand,  the  computed  phase  angle  between 
pitch  and  plunge  motions  and  the  flutter  frequency  were  not  significantly  affected  by  the  porosity 
parameter. 

The  main  conclusion  for  the  wind-tunnel  interference  problem  is  that  the  limit-cycle  flutter 
amplitudes  can  be  quite  sensitive  to  the  chosen  wind-tunnel  wall  porosity.  In  fact,  flutter  may 
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be  suppressed  completely  for  a  sufficiently  small  value  of  the  porosity  parameter.  Furthermore, 


the  free-flight  flutter  behavior  may  differ  substantially  from  the  behavior  found  in  a  porous  wind 
tunnel,  depending  on  the  chosen  porosity  and  blockage  ratio.  Consequently,  further  investigation  is 
necessary  to  assess  the  modeling  of  the  porous  wall  boundary  condition  and  the  correlation  between 
wind-tunnel  tests  and  free-flight  conditions. 
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APPENDIX  A:  DEFORMING  GRID  EQUATIONS 


Navier-Stokes  equations  in  non-dimensional  form: 


dQ  dF  dG  (dF  dG\ 

1k+!h+!k-Re  \lh  +  !k) 


(A.l) 


where  Q  =  (p,p«,p w,e)T . 

For  the  LHS  of  Eq.  (A.l): 


Q 


Q{x,z,t)  =  <2(£,£,t)  => 


dQ  _  dt,  dQ  ,  9£  dQ  ,  9t  dQ  _ 
FF~'3Fd^~l~'3Fd^~l~'drSz~ 


^tQ-^  +  +  Qx 


F  =  F{x,z,i)=F{lnta,i)=> 
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+  =  ^Gl  +  ^Gr 

dz  dQ  dz  dT  S  S 


Similarly: 
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dz 


&S  +  &t; 


(A.3) 


Substituting  Eqs.  (A.2)  and  (A.3)  into  Eq.  (A.l): 


£ tQ £,  +  CfGt;  +  £?x  +  ^ 


1  +  £z>%  +  Czty 


(A.4) 


Let  2  be  rescaled  by  the  Jacobian  by  defining  Q  =  J  lQ.  Then: 
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(J~lQ) x  =  J~lQ,  +  {J~%Q  ^Qt  =  -  ( J~l)zQ } 

(J~%Qk  =  =>  ^  -  (Z"1^)  5G] 

(Z-^G)?  =  J~%Q^  +  (Z^Cd  ?e  =►  ^  =  J[{J~%Q\-  (J~%kQ] 

(J~%Fk  =  +  (J~%kF  =>  ^  =  /[(/- -  (y-1^)^] 

=  r%F^  +  (z-*y  ^  =>  ^  =  J[{J-%XF\-  (. r%\F \ 
=  J~%G^  +  [J~%\G  =>  ^  =  J[{J~%G\  -  (. J~%\G \ 
=  J~%ZG^  +  {J-%Z\G  =>  CzG?  =  /[(Z-^G)^  -  (Z 
(z-1^  =  z-1^  +  =>  =  J[{J-%R\  -  (. r%\R ] 

(z- 1  y?)  ?  =  z- 1  y?c  +  (z~ 1 y  ^  =►  C**?  =  Z[(Z“ 1  y?)  ?  -  (Z- 1 1  y  c/q 

(Z_1^z5)^  =  Z-1^  +  =>  ^  =  J[(J~%S\  -  (J-%)^S] 

(J-%zSk  =  z-1^  +  (J~KzkS  =>  ^  =  J[(J-%fy-  (J-%\s\ 

Substituting  Eq.  (A.5)  into  Eq.  (A.4)  and  dividing  by  Z: 


(Z-'G)x -  (J-])xQ  +  ~  ( J-%\Q  + 

+  {J~%F\  -  {J~%\F  +  (J~ 1  C,XF) ^  -  {J~%\F 
+  {J~%G\  -  {J~%\G  +  {J-%G\  -  {J~%\G 
=  Re-{[{J~%R\  -  {J-%\R+{J-l^xR\  -  {J-X^X\R 

Rearranging  terms  in  Eq.  (A.6)  and  using  the  definitions  of  Q,  |  and 


(G)x-G[(Z  ‘)x+  (lf)^+  (Cr)^]  +  (IfG  +  lxZ1  +  fz<Z)^+(^Q+yr  +  yZ)t; 

-F[(y^+(tv)d-G[(y^+(yd 

^e_1{(lx^+|z5)^+ (y?+ y>)^  -R[(y^+  (y  d  _s[(yy  (y  dl 


(A.5) 


(A.6) 


(A.7) 
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Expanding  some  terms: 


(• J  )x  =  (x^  -  x^)T  =  x^zt;  +  x^t  ~  x^  -  x^x 


(A.  8) 


(If  )$  =  (x^Zx  -  XTZ^)^  =  X^Zx  +  X^z  ~  X^  ~  XzZQ 


(If) £  =  (xxZg,  -  X^Zz) i;  =  X^  +  XxZ^  -  X^Zx  -  X^Z^x 


Combining  Eqs.  (A.8),  (A.9),  and  (A.  10): 


(A.9) 


(A.  10) 


(•Ox+dfk  +  df^O  (A.  11) 

Important  to  note  that  Eq.  (A.  1 1)  is  known  as  the  Geometric  Conservation  Law  (GCL).  It 
will  be  derived  again  in  Appendix  B  to  show  that  it  actually  represents  the  conservation  of  volume 
as  the  grid  is  changed  in  time. 

Although  the  RHS  of  Eq.  (A.  1 1)  is  mathematically  zero,  it  is  not  the  same  numerically. 
The  assumption  of  orthogonality  is  no  longer  valid  for  a  deforming  grid.  Consequently,  one  should 
expect  some  error  in  the  solution  when  the  deformation  of  the  grid  is  relatively  large. 

Expanding  more  terms: 


(ick  =  (%k  = 
(Lk  =  (-Z5)c  = 


=*•  (I x)^  +  (Xx)t;  —  0 


(A.  12) 


(lz)^  -  ~  ~XQ 

(Cz)c  =  (^)c  =  -*K 
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=>&)!; +&)?  =  <> 


(A.  13) 


Plugging  Eqs.  (A.12)  and  (A.13)  into  Eq.  (A.7): 


(2)  x  +  (I  tQ  + 1 xF  +  +  (^2  +  txF  + 

=  Re~l[(i,xR  + 1 ZS)^  +  (C,XR  + 


(A.  14) 


Defining: 


F  =iQ  + If  +  id  R  =  lR  +  is 
G  =  t,Q+tF  +  i7G  S  =  +  l7S 

Substituting  Eqs.  (A.  15)  into  Eq.  (A.  14): 


(A.  15) 


Gx  +  ^  +  <%  =  Re~l{R^  +  Stf  (A.  16) 

Equation  (A.  16)  represents  the  Navier-Stokes  equations  in  non-dimensional  form  and  writ¬ 
ten  in  terms  of  the  computational  domain  variables.  It  is  important  to  note  that  the  Navier-Stokes 
equations  arc  valid  for  a  changing  grid  because  the  Jacobian  was  considered  to  be  a  function  of 
time. 
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APPENDIX  B:  GEOMETRIC  CONSERVATION  LAW 


The  velocity  of  the  grid  points  on  a  surface  S  can  be  described  by  W 5  =  (xx,zx)  f°r  ix-  z)  £  S 
and  the  velocity  of  all  grid  points  is  given  by  the  field  W  =  (xx,zT). 

Geometric  Conservation  Law  (GCL): 

^  f  dV=  f  Ws-dS  (B.l) 

dt Jy  Js 

Applying  the  Divergence  Theorem: 


-  f 

dt  Jv 


(B.l) 


But: 


dV  =  dxdz{\)  =7_1rf^(l)  =►  -f-  [  J~xd\d^=  /(V-W)  d\d^  (B.3) 

dt  in  in 


Expanding  the  term  (V  •  W) : 


i_1(V- W)  =  i_1 


r)iv 
3x  dz 


=  J~l 


/dt,du  dl^du  dt,  dw  d(^dw\ 

[dld^+d. ^+fcd^+dzKj 


(B.4) 


Rearranging  the  terms: 


i_1(V-W)  =i_1(V^-W^  + V^-W^) 


(B.5) 


But: 
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=>J-1V$-W$=  (/“1V^-W)^-(7-1V^)^W 

(B.6) 

(/_1vC-w)?  =  (y'vo^-w+y-^c-w^ 

=>  j~l vt;  •  w?  =  {J~l vt;  ■  w)?  -  (/_1  vg?  •  w 

Substituting  Eq.  (B.6)  into  Eq.  (B.5): 


y_1(v-w)  =  (/-1v^w)5  +  (/-1v^w)c-w-[(/-1v^5  +  (/“1vq?]  (b.7) 

It  is  also  known  that: 


•  W  =  +  §k)  •  (xx\  +  zTk)  =  E,xxx  +  %zzx  =  —J (x^Zz  ~  xxz(]  = 

'  W  =  (§i  +  §k)  •  (xx\  +  zxk)  =  t,xxx  +  t,zzx  =  - J{xxzi ■-  -  x^zx)  =  -£f 
Plugging  Eq.  (B.8)  into  Eq.  (B.7): 


7_1(V  •  W)  = 

(B.9) 

-W  •  {[(/-^.v)  5  +  (T-1?,)?]*  +  [(^)  5  +  (J~Kz)  dkl 
Rearranging  terms,  performing  the  dot  product  and  using  |  =  J~lt,  and  t,  =  J~lC,: 

/-i(v.w)  =-(|f)^-(y?-^[(i)^+(L)d-^[(|z)^+(yd 

=  -&k  - -  (*s)d  -  *[-(*&  +  Wd  (B10) 
=-&k-&k 

Substituting  Eq.  (B.10)  into  Eq.  (B.3)  and  assuming  that  the  Jacobian  varies  continuously: 
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f^{r')dW  =  |jaT(/-1)  +  (i)^  +  (t)d^^  =  °  (B.ii) 


Hence: 


axGr1) +  (£)*+ &)?  =  <> 


(B.12) 
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APPENDIX  C:  JACOBIAN  MATRICES 


A.  Matrices  A,  B,  and  M 

A  |  A  A  |  A  A 

It  is  necessary  to  know  the  matrices  A+,  A~,  ,  B~ ,  and  M  in  order  to  apply  the  numerical 

A  |  A  A  |  /V 

method  discussed  in  Chapter  II.  Matrices  A+,  A~ ,  B+,  and  B~  are  obtained  by  assuming  flux 
splitting.  These  matrices  arise  because  matrices  A  and  B  can  be  written  in  terms  of  an  eigenvalue 
factorization: 


A  =  LaAaLa  1  and  B  =  LBABLB  1  (C.  1) 

where: 


- 1 

0 

0 

o 

iS 

0 

0 

1 - 

O 

0 

^2 

0 

0 

and  Ab  = 

0 

Ef 

0 

0 

0 

0 

A-3 

0 

0 

0 

Ef 

0 

1 - 

O 

0 

0 

K  . 

0 

0 

0 

Also,  the  diagonal  matrices  A^  and  A B  can  be  split  into  the  summation  of  two  parts: 


=  +  and  kf  =  Jj?+  +  Xf-  (C.3) 

where: 

+  A,,  +  |X,|  Xi—  |  A,,  | 

K  =  — 2 —  and  ^  = — 2 —  (C4) 

In  doing  this,  matrices  A  and  B  can  be  written  as: 

A  =  Ea(A+  +  A-)E;1  and  B  =  EB(A+  +  A^1  (C.5) 
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1. 


Matrix  M 


0 

0 

0 

0 

AS'  _ 

m  21 

«i(l/p)c 

«2(1/P)C 

0 

nm 

«2(l/p)c 

«3(1/P)C 

0 

mi 

>1142 

"743 

«4(l/p)c 

where: 


m2\  =  -ai(w/p)?-02(w/p)5  n?3i  =  -a2(w/p)^  -  a3(w/p)^ 
n?4i  =  a4  [-e/p2  +  (ir  +  w2)/p\^-al(u2/p)^-2a2{uw/p)^-ai,(w2/p)^ 
>1142  =  -a4(w/p)£-m2i  »?43  =  -a4(w/p)^-m3i 


(C.6) 


(C.7) 


and 


<*i=p(5!2  +  C|)  a,  =  ^4  a3=,,(g  +  ig)  a4  =  ^  (g  +  g)  (C.8) 
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